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I. INTRODUCTION 


A. GENERAL PRESENTATION OF THE PHENOMENA 


Great progress has been made in the static and dynamic 
analysis of complex structures through the continued devel- 
opment of discrete element methods of structural analysis. 
Tremendous inprovenments in computing power have mad? 
possible the study of fully nonlinsar problems. 

The analysis of the rasponse of a structure submerged in 
a fluid, is severely complicated by the intrusion of signif- 
icant fluid-structure interaction 2ffects. Recently, the 
development of a variety of surfaca interaction approxina- 
tions has provided the means for a more efficient analysis 
of the coupling between the structure and the surrounding 
uid. 

Computer codes for structural analysis are well- 
developed so that the fluid-structure interaction is, for 
tha most part, handled by adding new capabilities to 
existing structural analysis programs. 

It is awell known fact that the primary threat toa 
Submerced structure is the shock wave that results from an 
underwater explosion. However, th2 complexity of <he phys- 
ical phenomena involved, along with the difficulties encoun- 
tered in obtaining experimental results have been serious 
drawbacks for the analysis of thess types of problems. But 
there is a definite need for investigations of large defor- 
mations and buckling problems in a structure submitted *o an 
underwacter explosion. 
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Be. OBJECTIVES 


This study deals with the nonlinear response of a 
submerged cylindrical shell +o a shock wave. ne sexes 2nd 
finite element code EPSA (Elasto-Plastic Shell Analysis) 
{Ref. 1] which includes nonlinear effects anda surface 
interaction approximation was selected for the study. In 
order to alleviate the tedious interpretation of results at 
points throughout the shell, PATRAN-G, a color graphics 
system, was used. PATRAN-S allows for a global representa- 
mem Of a quantity distribution over the structure rather 
than the discrete representation given by a computer output. 
The objectives of this study were to merge the finite 
element code EPSA with the color graphics syst2m PATRAN-G, 
and to conduct an analysis of che response of various types 
of cylindrical shells to a2 sphericai shock wave generated by 
an underwater explosion. 


a2 





II. THE EPSA COMPUTER PROGRAA 


A. PRESENTATION 


EPSA (Elasto-Plastic Shell Analysis) (Ref. 1] is a 
computer program deveioped by Waidlinger Associates and 
funded by DNA/NAVSEA/ONR for the purpose of the analysis of 
submerged stiffened shells under dynamic loadings. It incor- 
porates a number of specific features which are geared for a 
more efficient analysis. 


In particular, EPSA allows: 
“The analysis of shells in an acoustic medium, subjected to 
both low and high freguency shock loadings. 
“An efficient modelling of the elasto-plastic behavior 
-The inclusion of large displacement effects to analyze 
dynamic buckling situations and post-buckling behavior. 
-The modelling of stiffeners and internal structures. 
meme Cluid-structure interaction effect 

The following sections describe the equations of motion 
for a submerged structural shell in an infinite fluid 
sujected to a pressure loading and investigate the modelling 
Siecehe Surrounding fluid. The last sections are devoted to 
the finite element discretization as well as to specific 


features concerning EPSA. 


Be. EQUATIONS OF MOTION FOR THE SHELL 


Considering a thin shell of thickness h , volume V, area 
R submerged in an infinite fluid (figure 2.1). The shell 
stress resultants are defined from the stress tensor by : 


ee. 





h/2 ey 2 


Nis fase and ase | Oy4% Ae 
“h/2 sh/2 


Gime dasctribution of strains (2,5 °% 72 )>) is assumed <+o be 


linear, the curvatures at mid-height are (kjykr7k) . 


Od, (hee) TOP 


Diy (N22) Nii (Nee) 





= 7 
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Figure 2.1 Shell Stress State. 


Meeplying the principle of virtual work gives: 


rT i z 
[ ts¥teen aa -[ ©) {Su} dR - fs {u} {Su} dR (2.1) 
R R 


R 


Where 

{uj} = (U,,U5,W yr is the displacement vector at ¢ach point 
fep—(Ni1e-N22 -Ni2-M11 -M22,4%12 a is the stress resultant 
vector 

fe} = (11-25, 622, 4k), eee is the strain resultant 
vector 


Ois the mass per unit area for the shall. 


Th2 symbol § will refer =9 a virtual quantizy, and the dot 


or star denotes a differentiation with respect t> time. 
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miemenacst term of equation (2.1) represents the virtual work 
of internal forces , the second rapresents the virtual work 
of external fSrces (i.e. pressure, point loading, etc), and 
the third one represents the contribution of inertia forces 
in the virtual work. Thus, this last term expresses the 


effects of dynamic phenomena on th2 structure. 


C. FLUID MODELING 


In the case of a submerged structure, the external 
forces are the pressures applied at the fluid-structurs 
interface. 

As the shock wave hits the stuctuze it gets reflected, 
thus inducing a pressure term p, . In addition, the motions 
of the shell will also generate a radiated wave, Math a 
pressure contribution Phas 
Therefore, the pressure at the fluid structure interface is 


the sum of the incident, reflected and radiated pressures: 


P S Pi : Pr = Pra (2392) 
Where 
p = Total dynamic pressure. 
Pp; = Pressure associated with inclient free field pressure 
wave. 


Gee- Retiected pressure due to the iateraction of the inci- 
dent wave with the structure, the structure being fixed and 
magi d . 

Pra= Radiated pressure due to the normal movements of the 
surface of the structure in the fluid 

Pco= Prt+ Prais called the scattered pressure. 


Mae wwe-nods for getting the scattered pressure p,will 


now be investigated. Assuming a spherical wave in an 
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acoustic medium with a sound speed c, the pressure is 


determined by the well known wave 2qguation: 
2 Z 
aay Pp = 28 (2.3) 


and the proper conditions at the boundaries 2f the fluid 
domain. 

One alternative for solving th2 previous orcblem would 
Seeemmply to use a finite element discretization of part of 
the fluid domain, impoSing a radiation condition at the 
boundary [ Ref. 2]. However, this would reguire a very 
Sumi icant fraction of the computing effort that could not 
be devoted to the structural modelling. 

Therefore, a more efficiant way of computing the scattered 


pressure is required. 


The Doubly Asymptotic Approximation (DAA) imparts upon 
the structural model a surface loading composed of incident 
and scattered waves. 

In the high frequency limit, it can be shown that the 
Scattered nodal force vector {F,} is related to the wave 
particle velocities normal to the structure's surface by 
(eect. 3): 


“= CA J {U5} (2.4) 


Where [A] is the diagonal matrix of nodal areas (areas asso- 
Ciated with each node) and {Us} is the vector of nodal scat- 
tered normal velocities. Therefor2, in this high frequency 
case the shock wave is simply a plane wave and equation 
(2.4) states that the pressure is proportional to the fluid 


velocity. 
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inthe low frequency limit the fluid structure interac- 


tion is governed by the relation: 
* 
(F.} = (4,]{0,} (2-0) 


Where (u.” = S(0,} is the nodal normal acceleration vector 
amet mj] is the added mass matrix computed in an incompres- 
Sible fluid. 

Thus, in the low frequency cas2, the loading is due +9 
Seems wmotion of the frigi Structure in an incompressible 
fluid, a problem typicaily found in hydrodynamics. 


When a broader range of frequencies is considered, one 
can combine the two previous E€quatidns with the differential 


equation governing the scattered pressure (Ref. 3], giving: 


CAM (3 + oct, }# {Fg} = pele} (2.6) 


- d 
Where {Ff} Bee {F. } 
Defining the vector cf nodal scattared pressures {ps} by: 
{ps} = CA} (FP,} 
we get: 


CM UJ{ps} + pcCAl]{ps} = oc{ My] {Us} (27) 


Which is the set of equations that gives the scattered pres- 
Sure at 2ach node of the wetted surface of the shell. 

Equation (2.7) gives exact results in both the high and low 
frequency limits. Thus, DAA allows the modelling of the 
fluid-structure interaction via a coupled set of differen- 
tial equations at the wetted nodes of the structure instead 
of modelling the whole fluid with a finite element grid. 
The use of DAA will free some memory space in the computer 
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moree better modelling of structural behavior while az <zhe 
same time giving some reasonably accurate results for the 
loading of the structure. 

Therefore, the equations for the study of underwater 
shocks will consist of a coupled set of structural equations 
that come from the Principle of Vir<ual Work and of fluid 


equations at the wetted nodes that come from DAA equations. 


D. FINITE ELEMENT PROCEDURE 


:. iscretization 


imenmprtnc.zpue Of Virtual work is rewritten for the 


structure introduced in section B [Ref. 4] 


[er {$<} dk -f (p} {Su} dR +f ocaFtsuyae =B0 (2.8) 
R R 


R 


The surface of thé region is covered by 2 quadrilat- 
eral mesh cf N elements, each having an area A, . The 
previous integral can then be replaced by 2@ summation of 


integrals over A; .: 


N 
[ef {seyIR = yy (s}, {Se}. dR (2.9) 
i=1 A; 


The integrations over A. are then pexzformed by dividing the 
: z k 
element into four regions A; (figure 2.2) We have then: 


4 
k 
| {s}, {Se}, dR = 3 (s}, {d2,}; Aj (22010) 
A k=1 


3 
and therefore equation (2.9) becomes ; 


[st {Se3dR = \ i (s}; (oe, }. a. (2.11) 


R 1=1 k=1 


Using the sam2 procedure, it can also be derived 
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4 
it k G 
[roviss dR= ye {p} {5 uj A; = {(P} {Sq} (2.12) 
R ec 1 
. N 4 
om - 
[ecu d R= \ ) (ak {Su}; i = C4] {q} (9q} (2.13) 
i=] k=l 
Where [M] is the mass matrix , {P} is the vector of exter- 


nally applied forces, and {q} is the vector of nodal normal 
displacements for the structure. 


7 
| 


eG ee co cee SD ce ee a Ee ee 





—— 


Figure 2.2 Grid Points in EPSA. 


By definiticn, finite elament discretization can 
express the displacement {u} at any point of an element as a 
function of the displacements at the corner points of the 
element, defined by {fg}. 


{uj} = (HJ {a} (2. 14) 


Where [H] is 2 6x+t2 matrix of interpolation functions. 


Combining the derivatives of {u} will give the strain vector 


{oe}. In matrix form, equation (2.14) gives 
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ff) = ([8B] {q} (2.15) 


Where {B] is a 6,12 matrix function of the geometry of the 
element as well as a function of the previous deformaticns 
in the case of large displacements. 

Using the previous result, equation (2.9) is rewritten in 
the following way ;: 


T rT rz 
[o {Se} dR= \ , (s}, (8, J, (593, AL= {FY (6q (2.16) 
: ee eel 


Where {F} is the internal force vector. 


Combining the previous equations in equation (2. 8) oe ee 


principle of virtual work becomes: 


[M4] fq} = : (ee teh (2.17) 


1=1 


Therefore, the principle of virtual work has beén tranformed 
into a system of ordinary differential equations which ars 


much more amenable to numerical treatment. 


bn Ee SA each arbitrarily shaped quadrilateral 
element is defined by four corner nodes, each having three 
translational and no rotational degree of freedcm. In order 
to represent the behavior in bending Sight nodes not contig- 
uous with the element are also used (figure 2.3). Every 
element accesses twelve nodes and has twenty deqrees of 
freedon: three translaticnal degrees at the corner nodes 
and one corresponding to the displacement normal to the 
surface for each of the 2ight exterior nodes. The bending 
behavior (second derivative term) 1s expressed in terms of 
the nodal displacements via a finite difference technique 
(Ref. 4}. 
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Then the nodal displacement vector of element i is 


simply: 


{q} = Beano tage eo S84 eh ye oe MA pac ern) Y (2.18) 


Conventional finite elemant codes utilize +thres 
translational and two rotational degrees of freedom at each 
nod= (each element has 5 4=20 d.o.f. aS in EPSA ). However 
the masses associated with rotational degrees of freedom are 
very small, leading to numerical instability. The use of 
the aforementioned formulation alleviates this problem since 
only translations are considered and, im ead cien- the 
number of unknowns is reduced, leading to simpler and more 
efficient computations. 

It must be pointed out that in order to model the 
bending behavior at the edges of th2 shell, a set of artifi- 
Cial nodes has to be created. The finite element grid will 
then consist of the nodes iefined in the input file plus the 
artificial nodes along the boundary of the sheet 
@Migure 2.3). 


2. Strain Displacement Relation 


ThewepEinmcaple cf virtual work deals only with 
strains. Since the finite element approximation gives the 
displacement at each point, tha equations relating the 
strains to the displacements are nesded. In this study, the 
Donnell-Vlasov nonlinear kinematic 2quations of shell “theory 
are employed, and the strain-displacement relations are 
described in table I. 

Using equation (2.15) in the equations detailed in 
the previous section will give the finite element approxina- 
tion of strains in terms of nodal displacements. 


ZN 





TABLE I 


du dh 





2k * ae tee Ed 
eZ h, 96) ho %S> hh 2 35. 


aw 
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Where =x = 
o5 
ae are local coordinates 


h, -ho are scaling coefficients. 
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3. Shell Constitutive Relations 


The shell constitutive relations relate 


Tesultant rate vector to the shell strain 


Matrix terms: 


Ze 


-Donnell-Vlasov Shell Equations 
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{fs} = (D]@} (2. 19) 


Where [D] is the Elasto-Plastic tangent stiffness matrix. 


ay 





Se 


Fagure 2.3 Nodal Points Organization. 


The stress-strain relation used in EPSA differs fron 
the classical Elasto-Plastic theory in that the formulation 
involves shell stress resultants rather than stresses at 
points throughout the thickness of the shell. In other 
words, EPSA uses relation 2.19 integrated over the thickness 
of the shell. Thus, there is no need to compute the stresses 
throughout the shell, which results in a significant savings 
in storage space and processing ¢+ims. However, the stress 
resultants Nj and Ms 4 cf the shell theory are not sufficient 
to describe the state of stress, and certain higher order 
moments mus*= be combined with the stress resultants to form 
the dynamic variables of the problem. The coefficients for 
the integrated constitutive equation have been determined 


uSing experimental results [Ref. 5]. 
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Meee COnStitutive equations consist of a yield sondition, 3 
strain hard2ning law anda flow rule: 

-The yield condition indicates whether part of the shell has 
Started to yield. (figure 2.4) 

-The strain hardening law gives th2 evolution of stresses in 
the shell after plasticity is reached. 

-The flow rule gives the plastic strain rate in the shell 
arcer plasticity. 


ee ee eee 


Socom 





Figure 2.4 Yield Situation in the Shell. 


4. Solution Procedure 





EPSA uses an explicit central difference scheme to 
solve the virtual work and fluid loading equations detailed 
in section B. As indicated in appendix B, the explicit time 
integration procedure requires a small time step and is only 
Senditionally stable. However when stable, it always 
converges to the exact solution, as opposed to implicit 
schemes that are unconditionally stable but may lead to 
unrealistic results. Furthermore,in problems involving the 
treatment of shocks, accuracy requiremenzs preclude the use 
of large time steps. The central difference scheme seems, 


therefore, particularly optimal. 
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Assuming we know the solution at tine = ,the central differ- 
ence scheme applied to equation (2.17) gives th solution at 
time t+At : 


tet t z 
= «+ At CEPR ~ 2 Ph) (2. 20) 


Mi 


Where M, is the mass of node i, {P}; re the externally 
appied forces and (F}, are summed over all the elements k 


foaming nods i. 


fae tormmlaticn of the equations is in the initial 
configuration and all equations are solved in the initial 
geometry in accordance with the total Lagrangian formulation 
{Ref. 6]. 


E. EPSA CAPABILITIES 


The structure to be analyzed is conceptually divided 
into constitutive parts called "sh2ets." Each sheet is a 
curved section of a shell with an arbitrary number of nodes 
and elements (figure 2.5) Its shap= is limited to a surface 
that can be described by a smooth continuous function 
without discontinuities in its slope. The elements within 
the sheet are limited to a rectangular organisation (figure 
ee) - 

Thus, e@cylinder with end caps would consist of *hree 
Sheets: a cylindrical sheet anda sheet for each end 
Meegure 2.5). Three sheets are required because of the 
slove discontinuity at the edge between the cylinder and the 
end caps. 

Dividing the structure into sheets isolates the diffi- 
culties associated with the boundaries into a set of artifi- 
cial nodes along the perimeter of the sheet. When several 
Sheets are required ,EPSA takes car2 of the compatibilities 


of displacements, rotations and moments along the edges. 
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Figure 2.5 Sheet Organization. 


Thus, any arbitrarily-shaped structure can be analyzed 
using EPSA by dividing it into a number of sheets. 

Two options for choosing elements are available in EPSA. 
The first option exists to employ 2 generalized quadrilat- 
eral element. The second option exists to employ a rectan- 
gular element and uses a specialized formulation for this 
type of element. 


The coordinates which can be a2ither cartesian or cylin- 


tQ 


drical always lie within the plane of the sheet. The 
direction lies in a positively outward direction normal to 
the sheet. Each sheet contains its -own local coordinates 
system, there is no global coordinate system when multi 
sheets are merged (figure 2.6). 

Pricer to generating a finite 2lement mesh for a sheet 
one must establish the sid2 numbers of the sheet. The side 
numbering scheme is arbitrary as to the choice of sheet 
humber one. Hcwever the specification of sides 1 +o 4 nust 


proceed in a counterclockwise direction when the sheet is 
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Figure 2.6 Coordinate System. 

viewed from the positive z direction. At the intersection 
of side 1 and side 4, element 1 is first defined. Then the 


node/felement number is incremented by one until it reaches 
side 2. Then it returns to side 4 and continues the count 
for the next row of nodes/felements. 

Thanks to the exclusive use of quadrilateral elements 
and to the specific numbering procedure, the table of 
connectivity is implicitly defined when generating the nodal 
points, thus simplifying considerably the input require- 
ments. 

The inclusion of structures internal to the cylindrical 
Sheet is enacted in EPSA through the use of internal sheets. 
Structures internal toa cylinder are therefore nodelled as 
individual sheets having the sane capabilities as any 
general EPSA shest. 

Two types of internal structures are available: 

-Sheets (beams, plates or shelis) whose connection to the 


cylindrical shell lies parallel to the axis of the cylinder. 
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-Sheets (beams, plates, shells) whose connection lies in the 
circumferential direction of the cylinder. 

Furthermore, in crder to moiel internal equipmen+ 
(machinery, etc ) Which does not ieform but contributes to 
the inertia of the system, concentrated masses can be input 
at nodal points. 

The user aust be aware that the previously discussed DAA 
is only implemented for cylindrical stuctures. Prior «o the 
finite element analysis the user must compute the added mass 
(virtual mass) matrix defined in 2quation (2.4). Theses 
done by using *h2 ACCESION program, which creates a virtual 
mass file that EPSA reads when computing the fluid-structure 
interaction. 

Finally, EPSA in its latest version takes local cavita- 
tion into account. When the total pressure computed by EPSA 
is found to be negative, j2eisomesd IDLY Set 25> —2er0 Since 


water cannot withstand any tension. 


F. USING EPSA 


The input file for EPSA can be generated either 
directly, or via the interactive program PREPSA that prompts 
the user to give the input data. When the input file is 
created, all the data are in free format. 

The nodal points can be generated semi-automatically 
(see the user's manual), and this option is very helpful and 
time saving when generating big models. EPSA can be run 
interactively for small models or on batch for bigger jobs. 
For instance, a cylindrical shell with 1440 elemenzs and 
1517 nodal points takes 0.0129 sec. of VAX CPU per time step 
per element. The whole model would take about half an hour 
for 100 steps. 

EPSA creates an output file in which all input data is 
echoed, and outputs the pressures, stresses, strains, veloc- 
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ities, displacements at the nodal points requested by the 
user in the input deck. 

The value of the time increment At should be selected so 
jmat : 


me < 1/26, ( 07z)° Cee 
Whereo_ . is the smallest distance between th2 midpoints of 
Opposite sides of an element, for all elements of the sheet 
(figure 2.5), and 1/2 is a safety factor. Etmozher words, 
the time step increment has to be Less than the time taken 
by a wave to propagate from an element to another. In equa- 
tion 2.21, (E/o) is simply the wave speed in the material. 

The virtual mass array (VMA) is created on unit 20, 
therefore one should not use this unit for any other purpose 
than READ operations. 

Because of the simple organization of its input file, 
EPSA has been found fairly easy to use. Thé user can perform 
Major changes in the model very quickly and effirtiently. The 
ACESION modul2 has been found to work well except for cylin- 
ders of large dimensions (L=1400", D=240") for which the 
size of +he virtual mass array grows unexpectedely froma 
reasonable 4 blocks to 190 blocks of VAX memory size. 

However, EPSA has been design2i for some specific types 
of fluid-stucture interaction analysis and its limitations 
should be pointed out: 

- Only beam type stiffeners can be considered 
= fhe fluid structure interaction is only enacted for 3 
Circular cylindrical sheet, which takes away much of the 


flexibility the progran. 
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III. BPSA/ PATRAN-G INTERFACE 


A. INTRODUCTION TO COLOR GRAPHICS SYSTEMS 


In dealing with the respons= of a structure to 


jw 


loading, quantities such as stress2s, s*rains, velocities 
and displacements must be analyz2i ata number of nodal 
points, which makes the interpretation of computer outputs 
very tedious. Color graphics systems offer an effective 
solution to this problem by providing a global representa- 
meen Or a quantity distribution over a stucture rather than 
a discrete representation given by a comput2r output. A 
color graphics system consists of an interface between the 
computer, the us2r and the color terminal. A typical system 
would be a software package that allows the user to creat? 
models on the screen as well as to display any data in terms 
@=e color intensity. Deis Known that avpicture is worth 
several hundreds of words. Thersfore, merging a finite 
element program with a color graphics system would be a 
major improvement in engineering analysis. 


PRURAN-G (REE. 7] is a color graphics system specifi- 
cally designed for finite elements, it permits the engineer 
and the computer to work together towards the creation of a 
mod2l. The designer creates an image on the screen and 
PATRAN automatically translates th2 physical model intoa 
standard finite element input deck. Another important 
feature of PATRAN-G is its ability to assist the user in the 
interpratation of results through its post-processing facil- 
ities which include: the deformed geometry with magnified 
deformations,the color coding of elements based upon any 
response parameters such as displacements, stresses and 


strains. In particular, the contour levels of the 
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aforementioned quantities can be Superimposed on 2 
3-dimensional image of the model, thus allowing for a global 


analysis of a complex structure. 


Be. MERGING EPSA AND PATRAN<-G 


As described in chapter I , ‘the structures under study 
have fairly short EPSA input files. FUGthOEmoOne, iin this 
study dealing mith  fluid=structure interactions on 3 
circular cylinder, only structures that consist of one sheet 
are considered. Miereroce, —~ because of the simplicity of 
both the input file and the stuctur= under study, there was 
no need to design a module converting a PATRAN-G model that 
ils created on the screen into an EPSA input fils. 


The remaining tasks were the following: 


-Display the original finite element model defined in the 


EPSA input file on the screen (original geometry). 


~Display the nodal points and :23lement results that are 


computed from EPSA on the screen (postprocessing). 


The input of a finite elanent model into PATRAN-G 
can be done by creating 2 "neutral" file, as described in 
PATRAN-G terninology. The neutral filet is intended to 
provide a simple link between PATRAN-G and the outside 
mOoLld. It is written entirely in 80 character card images 
gag all the data is organized in snall "packets" of two or 
more card images. Bach packet contains the data for 43 
fundamental unit of the model such as node coordinates or 
elements definition. Since our only purpose was tc display 
the original geometry of the structure, a limited number of 


lAdditional information about nautral files can be found 
in the PATRAN-G user's manual (Ref. 7] 


31 





data packets (4) was to be created: 


-File <itle (packet 25) fwewerca tas COhtalhing the user 
title. 


“Summary data (packet 26) : two cards containing the number 
of nodes, elements andthe date andtime at which the 


neutral file was created. 


“Node data (packet 1): contains all information concerning 
nodes: node number and coordinates in a global coordinate 


frame. 


-Element data (packet 2) : contains the connectivity table 


for the finite element model. 
meena, OL file (packet 99) : end-of-file cards. 


We have seen in the presentation of EPSA in chapter I «hat 
the nodes are defined in a local, sheet-attached ccordinatse 
system, that artificial nodes are craated along the edges of 
the sheet to nodel the bending behavior, and that nce connec- 
tivity table was input. Therefore, the translator module 
that would be created had to: 

-skip the set of artificial nodes and properly renumber the 
ra a 

-perform a change of coordinates for all local data 
-generate the connectivity table. 

It was decided to employ a modular design in which 
each routine would perform a specific task. A modular design 
would allow further changes to be mad? quickly and ¢ffi- 
Clently by modifying only the relevant routines. The imple- 
Mentation in EPSA was made using a series of "calls", «hus 
Minimizing the risk of interference with the finite element 
computations. 

The translater module cr3ated was made of four 


routines: 
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~PRELIM >: computes the number of nodal points, the number 
of elements and displays the first two data packets (25,26) 


aot OTF >: scans through the aod2s, skips artificial 
nodes, renumbers the grid, performs the required changes of 


coordinates and dispiays the node data packet (91) 


-SHCONN >: scans through the nodes, connecting each node 
to the elements it belongs to and displays the element data 
packet (02) on the neutral file. 


-ENDNEU : displays the end-of-file packet (99) 


The translator calls were inplemented in the routine 
REPORT of EPSA. Any run of EPSA creates a neutral file on 
unit 19. The neutral file name is therefore FORO19.DAT if no 
"ASSIGN" statement has b2en issuei pricr tothe computer 
run. The finite element model might then be displayed on the 
graphic terminal (Ramtek, Tektronix) via the neutral input 
mode of EPSA (see {[Ref. 7] for more details). An example of 
finite element model output on Tektronix 4014 is given on 


Eeoure 3.1. 


3. Implemantaticn of Postprocessing Capabilities 


Postprocessing capabilities exist to assist the user 
in the interpretation of computer results. It is simply 3 
process of generating displays and reports based upon a 
combination of geometry and the results of an analySis. 

The results cf analyses ar2 transmitted to PATRAN-G 
for postprocessing via "neutral rasults files" as decribed 
in PATRAN-G terminology. Unlike the model neutral file 
described in the previous section, results files are in 
binary rather than in card image forn. 
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Pigure 3.1 Example of Finite Element Model Display. 


One can distinguish between two major kinds of post- 
processing displays: deformed geometri2s and elsment 


quantities. 
a. Deformed Geometry 


A displacement results data file required by 
PATRAN-G had to be created. Again, the module created had to 
Skip the artificial nodes, perform changes of coordinates as 
well as to write the nodal diplacements in the neutral 
results file. The displacement results data file was 
created in module NEUDISP. Its organization 1S given on 
maple II. 

A small modul2 PLOTDISP that decides at which 
time steps the results should be output was created. The 
"Call" to PLOTDISP was implement2?i in module COMPUTE of 
EPSA. The overall structure of the translator module is 
presented on table IV at the end of the chapter. 
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TABLE II 

| ‘ : : 
| Organisation of Displacement Results File 
| 


rigid body mode. 
=. Velocity normdl to the shell 


| 
| 
Column Content | 
| | 
| ! 
{ 1 X displacement in global coordinates | 
| 2 Y displacement in global coordinates 
| 3 j "a 8 " 
4 Displacement normal to the shell without 
| 
| 
| 
| 
[Ss 


a A i, Ee a ES 


A sample of deformed geometry output on Tektronix 4014 is 


Meven in figure 3.2. 
b. Element Quantities Postprocessing 


Any element related quantity can be the target 
of postprocessing. In general thes2 types of quantities are 
the results of finite element computations supplied to 
PATRAN-G through the neutral element results file. The 
neutral element results file is different from the neutral 
displacements results file detailed in the previous section, 
however it shares a Similar format "Ref. 7]. Element quan- 
tities such as stresses in local and global coordinates and 
von Mises stresses are computed in module NEUSTRE whose 


overall structure is similar to NEUDISP described earlier. 
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Figure 3.2 Deformed Geometry Output. 


The organisation of. the neutral alement results file is 
given on table III. 

As described in chapter I, EPSA does not ccmpute 
the stresses through the thickness 9f the shell. Instead one 
uses the integrated quantities of the shell theory [Ref. 8] 


h/2 Me 
vase [: dt and nse [3 tdt 
=n /2 


One can expect the stresses on the shell to be 
Maximum at the extreme fibers. NEUSTRE computes the von 
Mises stresses at the top and bottom fibers and writes the 
Maximum value in the neutral file. At the surface of the 
Shell no shear stresses are involved. Assuming a linear 
distribution of bending stresses and a uniforn distibution 


of membrane stresses, one can easily derive : 
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Ci5 = My /h + Ms 6/h?e (a) 


The first term of the previous equation is the membrane 
force contribution, the second is the bending moment contri- 
bution. The von Mises stresses are then computed using the 
well-known relation: 


Oo = ( Oy mor 0, + 0, ) (3752) 


In a Similar way ‘+o the displacements results 
Eve, a module PLOTSTRE that desides whether or not +9 
output the element results was created and called from 
COMPUTE. The Sverall structure of the translator is given on 
table IV. 


| 
: 


TABLE IIl 
Organisation of the Neutral Element Results File 





| 
| 
{ 
| 
Column Content Description | 
{ 
Z2 stre,1 Element local stress, direction i | 
23 st re ' 2 tf 19 " 1f 5 i 
25 stre,4 Element global stress, direction x | 
26 stre,5 u od " o ~ Yau 
{ 
27 stre,6 Element global shear, direction xy | 
31 von von Mises stress . 


—E—— 
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CoeUVopeayitugeene Results Quantities on the Screen. 


The EPSA input deck has been modified so as to 
create the results files (displacemants, elements) required 
by PATRAN-G. At the end 9£f the second input card the user 
specifies the number of displacement results and the number 
of element results files to be created (at least one). 
Obviously, those two inputs are also in free format. The 
neutral results files (elements, displacements) will be 
generated at equal time intervals as ths computation 
proceeds. The results corresponding «to the last time step 
are always output. 

The element results file is created on unit 16 
and the displacement results file on unit 18, thus corre- 
sponding to files FORO16.DAT and FOIRO18.DAT raspectively. A 
hew version of those files is created each time an output is 
requested. 


Example: 


If five (5) aeutral element results files are requested on 
the input card of a run of 20 steps, five files FORO16.DAT;1 
to FORO16.DAT;5 will be created, corresponding tS time steps 
4“ to 20 respectively. 

For the displaying of 2lament and nodal points 
results, the user will cefer to [Ref. 7] section 20. The 
mete Of the run (first card of EPSA input deck) will be 
displayed on “the screen along with the time at which the 
resuits have been reguested. 
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IV. DESCRIPTION OF THE UNDERWATER EXPLOSION 


A. PRESENTATION 


An explosion is a chemical or nuclear reaction in a 
substance (water) that converts an original material into 
other products plus a significant amount of 2nergy. The 
process of th? explosion produces vary high temperatures and 
pressures and occurs with extreme rapidity. As the result of 
the explosion, the initial mass of explosive becomes a very 
hot mass of gas at tremendous pressures; it is then obvious 
that these conditions will affect the surrounding medium. 

The fact that the water is compressible leads +o the 
conclusion that the pressure applied at some location in the 
liguid will propagate through it as a wave disturbance 
{Ref. 9]. It must be pointed sdsut that the pressures 
involved in an underwater explosion are usually so large 
that the wave velocity cannot be assumed independent of 
pressure. Thus, the small amplitud> wave theory detailed in 
{Ref. 9} does not apply and this will be the source of many 
complications in describing the behavior of the shock wave. 

The first cause cf disturbance to the water in an under- 
water explosion is the occurrence of the pressure step at 
the water boundary. Immediately upon its arrival, the pres- 
Sure begins t> be relieved by an intense pressure wave and 
Outward motion of the water. Por explosives like TINT or for 
huclear explosions, the pressure rise can be considered as a 
discontinuous step, andis then followed by a roughly expo- 
nential decay. The duration of the phenomenon is of the 
order of a faw milliseconds (more for nuclear explosions). 
Once initiated, the pressure disturbance is propagated radi- 
ally outward as a compression wave, also called a shock wave 


because of the steep pressure step at its front. 
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Ouome to the explosion, the velocity of the wave 1s several 
times the "“acoustic"™ speed of tha small amplitude theory 
(c=5000 ft/sec) ; but approaches this limit very rapidly as 
the wave advances outward. 

The pressure level in the wave falls more rapidly with the 
radial distance than what is predicted with the small ampli- 
tud@2 theory, but approaches this behavior at large 
distances. 


B. BUBBLE EFFECT 


As a result of the explosion, the initial mass of explo- 
Sives becomes a hot mass of gas 2t tremendous pressures. 
After the principal part of the shock wave has been emitted, 
the gas pressure is considerably decreased but is still nuch 
higher than the squilibrium pressur2. The water in the imne- 
diate region of the sphere or "bubble" of gas has a large 
outward velccity andthe diameter of the bubble increases 
rapidly. The expansion continues and the internal gas pres- 
sure decreases gradually, but the motion persists because of 
the inertia of the ocutward flowing water. When the gas 
pressure falls below the equilibrium value, the pressure 
defect brings the outward flow to a stop and the boundary of 
the bubble begins to contract at an increasing rate. fhe 
inward motion continues until the compressibility of the gas 
reverses the motion. Thus, the inertia of the water 
together with the elastic properties of the gas provide the 
necessary conditions for an oscillatory system. The oscilla- 
tions of the gas sphere may persist a number of cycles, ten 
or more oscillations having been detected in favorable 


cases. 
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Cae seer ACE EFFECTS 


In the case of explosions occurring at shallow depths, 
surface effects will complicate th2 aforementioned sequence 
of events. When the shock wave hits the surface,the atmos- 
phere cannot supply appreciable rasistance by compression. 
AS a result, a reflectsd wave with a negative pressure 
satisfying the zero-pressure condition at the surface is 
formed (figure 4.1). Thus, the resultant pressure observed 
is the sum of the direct and reflected pressures. Therefore, 
the reflected wave arriving at point A will create a sudden 
drop of the pressure to 3 smaller value. This is known as 
the "cut-off" phenomenon, typical of free surface effects 
(figure 4.2). 





| 
L 


Figure 4,1 Surface Effect on a Shock Wave. 


D. PRESSURE DETERMINATION 


As detailed in a previous section, the pressure decay at 


any point is roughly exponential so that it can be written: 


P(A) = By (A) e tf (4.1) 
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It has been found that the fundamental descriptors cf an 
underwater explosion attack are the charge size (equivalent 
weight of TNT) and the charge standoff (shortest distance 
between charge and target). Theoretical developments about 
the spherical wave detailed in [Ref. 10] have shown that the 
peak pressure Php at any point can b2 reasonably approximated 
by: 


Ry 
P, = kK, (W mR (4.2) 
where W is the charge size in pounds of TNT and R is the 


standoff distance in feet. 


It has been shown as well: 
Vea 5 A2 
6 = KW (W /R) (4.3) 
K, 0K, A, 7A, are empirically determined factors that depend 


on the type of explosives used. Their values for several 
types of explosives are given on tapdl2 V. 
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Figure 4.2 Cut-Off Phenomenon. 


4 3 





—— ED a i A Gy SEUSS es-ES gp) «SSS 


TABLE V 
Explosion Parameters 


cicinao lal 





7 


Ee. THE EXPLOSION IW EPSA 


EPSA features two different way of describing underwater 
explosions: 
“A discrete form in which the user inputs the pressure 
history at a finite number of times. 
“A functional form that uses equations 4.2 and 4.3  .Thea 
program requests then th2 various coefficients and parame- 


ters describing the explosion. 
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Pigure 4.3 Incident Pressure Decay. 
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The presence of free-surface effects can also be 
Meeouneed Lor With the input of a cut-off time after which 
the incident pressure is set to zero (figure 4.3). 
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A. MODELS STUDIED 


In order to compare the results of the numerical anal- 
ysis with the experimental data, all attention was focused 
on the Explosive Power Mester (EPM) model for which field 
tests had been conducted. The EPM model is a stiffened 
cylinder with end cars whose dimensions are given in figures 
ee 
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Figure 5.1 Explosive Power Meter. 
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Rather than modelling the end caps with additional 
sheets , it wae Geownd to be mors efficient +> taks into 
account the behavior of the end caps using two rigid end 
meocks (figure 5.1). The effects of the explosion on «he 
deformation of the end blocks is negligible compared «+o <¢he 
deformation of points located outside of the end blocks. 
Therefore, it was assumed that the aotions of the end biocks 
were pure rigid body displacements. 

In order to gain some insight into the influence exerted 
by the stiffeners and the end blocks, a preliminary analysis 
waS conducted on a ring stiffen2d cylinder without end 
blocks as well as on an unstiffened cylinder . 

In addition to the study of the EPM model, the influence 
that the location of the stiffeners had on the deformations 
throughout the shell was evaluated. By performing a compara- 
tive analysis of the deformations, it was intended to mini- 
mize and control the damage caused to the structure. 

The cylinders tested were subjected to an explosion 
occurring at the distance R= 200" from the cylinder. The 
location of the explosion was symmetrical with respect to 
the longitudinal and transverse axis of the cylinder (figure 
Bez) « A spherical type TNT explosion of intensity W=50 1b 
was selected for the study. It was therefore determined by 
the following parameters (chapter IV) : 

A, = 1.18 A> = -.185 
Ky 22505 Ko = .058 


The symmetry with respect to the xy and yz plane has 
been taken advantage of by modelling one quarter of thse 
model. After a certain amount of sensitivity analysis was 
performed on simple grids, a finite element grid consisting 
cf 1517 nodes and 1440 elaments was selected (figure 5.3) . 
For each of the cylinders studied, the time step chosen was 
equal to at= 3.106 s. The explosion process was studied 
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Figure 5.2 Explosion Location. 
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Pigure 5.3 FEM Discretization 1517 nodes, 14840 elements. 
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over a time interval of 800 time steps, that allowed for 2 
shock and after-shock analysis. 

The displacements computed by EP®SA consist of a combina- 
tion of rigid body displacements ani pure deformations. The 
deformation modes are of significant interest since they may 
induce buckling and even lead to the destruction of the 
structure. As mentioned 2arlier, the very stiff end blocks 
have pure rigid body displacements. The displacement of 
each node of the end block was subtracted from the displace- 
ment of each node cf the corresponding frow, giving the 
component of the displacement due +5 pure deformation. 

For each of the aforementioned cylinders, the deformed 
geometry and the coler-fillead contour plots of von Mises 
stresses as well as normal displacement were displayed. 
Using identical color assignments, the progressive aross 
evolution of the previous quantities were evaluated so as to 
allow for a comparative analysis of the evolution of phys- 
ical parame*ers throughout the shell. 

Color-filled contour plots allow for 2 global representation 
of a quantity distribution and have been found extremely 
valuable in the interpretation of the results. 

For printing and processing reasons, it was not possible 
to include color pictures in this document. Instead, the 
contour plots of the physical quantities under study have 
been included. 


B. ANALYSIS OF RING STIFFENED CYLINDERS WITH END BLOCKS 


1, ERM Model 


{ty 


The contour plots of von Mises stresses at time 
steps 20, 60, 100, 140, 180, 200 are provided on figure A.1 
momA.6 . As the shock wave hits the structure, i* appears 
as if the stresses propagate through the shell and reach 


their maximum value fairly quickly in 50 time steps. After 
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100 time steps, when the structure is fully enveloped by ths 
shock wave, the region close to the end block bécones 
heavily stressed (figure A.4). In addition, some concentza- 
tion of stresses at the locations of the stiffeners can be 
seen (figure A.5). At later time steps, th? pressure becomes 
decreased and there is a relaxation of stresses. However, 
*h2 region close to the end block as well as some spots 
located around the stiffeners remain heavily strassed, which 
may indicate local buckling (figur2 A.6). 


2. Controlled Deformations 


ieecrach ss to Obtain a more uniferk distributicen of 
displacements, the stiffeners have been shifted towards the 
end blocks. The time history avolution of the displacements 
was studied over an interval of 809 time steps for the EPM 
model as well as for the model with shifted stiffeners 
called EPM2. The ccntour plots of normal displacements at 
time steps 200, 400, 600, 800 for both models are provided 
on figure A.7 to A.14 

The EPM aodel shows a significant concentration of | 
deformations occuring, 2ven at late time steos (figure 
AetO), indicating a possibility of buckling. Although unex- 
pected, the fact that the region close to tha erd blocks 
undergoes the most severe deformations has been confirmed by 
experimental data. A possible explanation to this phenomenon 
is that the inertia of the cross-saction of the cylinder is 
relatively uniform along the cylinder, except at «he end 
blocks where it jumps to a much higher value. This disconti- 
nuity in the inertia results in concantrations of stresses 
that cause the aforementioned phenomenon. 

The cross-section inertia of the EPM2 model 
increases more smoothly because of the distribution of stif- 
feners along its axis. Thus, the concentration of stresses 


has a lower magnitude and the region clcse to the end block 
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suffers less damage than in the EPM case. It can also be 
seen that the deformations are more uniformly distributed 
along the axis of the cylinder. Above all, the deformations 
at the late time steps are not as large, indicating that the 
chances of buckling are significantly lower for the EPM2 
model (figure A.14). 

Therefore, by performing an optimization of the 
location of the stiffeners, the designer can counteract «he 
concentrations of stresses and tha buckling phenomena that 
occur in the régisn close to the end blocks. It is believed 
that controlled defcrmations can have a ver Significant 
influence on the survivability of a structure submizted «+o a 


shock wave. 


C. ANALYSIS OF UNSTIFFENED AND RING-STIFFENED CYLINDER 


It was decided to study the progressive gross responses 
of an unstiffened cylinder as well as that of a ring- 
stiffened cylinder without end-blocks. Both cylinders have 
the same external dimensions as the EPM model. Liew fang — 
stiffened cylinder is similar to the EPM model except for 
the fact that the end-block has been replaced by a standard 
stiffener. The evclution of von Mises stresses at time 
steps 40, 80, WOON tS prowided in Eigures A. 15 thrceugh A.17 
for the unstiffened cylinder. The evolucion of von Mises 
stresses at time steps 40, 80, 100, 150 is provided in 
figures A.18 through A.21 for the ring-stiffened cylinder 
without end-blocks. 

For the unstiffened cylinder itis observed that the 
stresses propagate quickly throughout the shell and that 
within a hundred time steps an instability phenomenon occurs 
showing the existence of local buckling (figure A.17). At 
later time steps the buckling spreais over the entire stiff- 
ened shell. 
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The evolution of von Mises stresses in figures -A.18 
mepeuwqn A.21 shows that the fring-stiffened cylinder can 
withstand much higher stress levels than the unstiffened 
shell. This result was expected since the stiffeners provide 
the stiffness required to withstand higher loads. At time 
step 100 the unstiffened shell is already subjected to local 
instability characterized by a “hard spot" in the middle of 
the model (figure A.17). On the other hand, the stresses in 
the stiffened cylinder are much more evenly distributed 
throughout the model, with high amdunts of stresses concen- 
trated around the locations of the stiffeners (figure A.21). 

It can also be observed that 4&4 significant concentration. 
of stresses occurs at the extremities of the stiffened shell 
fergure A.21). Recalling that the end-block has been 
replaced by a standard stiffener, the cross-section of the 
shell has a greater inertia at ths extremities due to the 
fact that the two stiffeners locatad at the extremities are 
close to each other. Therefore this phenomenon is similar to 
the one encountered when studying the EPM model. However 
tee Concentration of stresses for the ring-stiffened cylin- 
drical shell is less significant than for the EPM model, due 


moma smMaller discontinuity in cross-section inertia. 


Siz 





VI. CONCLUSION 

A FORTRAN module that merges the finite? elemert code 
EPSA with the colcr graphics system PATRAN-G has_ been 
designed and succesfully completed. The non-linear elasto- 
plastic responses of various types of submerged cylindrical 
shells have been evaluated using the EPSA/PATRAN-G systen. 

The analysis of the progressive gross responses of 4a 
ring-stiffened cylindrical shell with and-blocks is believed 
to provide useful information regarding the behavior of 4a 
submerced structure subjected to an underwater explosion. 
The influence of the location of the stiffeners on the 
deformations has been studied and is also believed to be of 
Significant help in the determination of an optimal design 
that will minimize the damage due to an underwater 
explosion. 

The utilization of the color graphics system in the 
interpretation of the results of analysis has been found to 
be an extremely valuable tool. It is the author's belief 
that the use 9f color graphics systems will become increas- 
ingly important in the analysis of complex phenomena such as 
underwater exploslons on submerged structures. 
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Von Mises stresses, time step 200, EPS model. 
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Figure A.9 
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Von Mises stresses, 


Figure A.15 
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Von Mises stresses, 


Figure A.16 





we ee SE coe OEP te Om ee OR ee SS ee | a ee ee Oe wae SS eae = eae ee eee Oe EE eg GE a ae 


piece 


N=" OVET ~~ 


£=' B26 


1='OOPTI ia 
H="OOTET ' Wis 


D='G08bI 
4=' 99591 TES 2) 

3 
3='@0Z81 Sa 


G='aaéét 


J='BO91Z i iaaecices ee 
El 


G='ABEEZ 


! 
Y=" 80052 


| ee Se ae SE Se SE gee SE age epee SE es SP ee ee SS age OSE gees SP gee ees SG 





oS Ge eer SSS ee CP SE SE gee CE ee ce Se eee ee eee ee pee SE ee 


6O-32966666 2 
TISHS WIITYONT TAD G3NSAITLSNN 


s 2 


ee 


é 


ee 


J 


isa 
a 


Ce eee ce en crn mc mam eae cence atc ee ec cea NN DE ae 


SS a Se ee EE ee ee EE ee EO i 


unstif. 


shell. 


Step 100, 


Yon Mises stresses, 
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Figure A.21 





_- ABPENDIX B 
REVIEW OF NONLINEAR FINITE ELEMENTS 


A. INTRODUCTION 


This appendix is irtended to give the reader some 
insights intd nonlinear finite elenents. Tne reader is 
assumed to have some previous knowledge of finite slomenz 
theory. The basic principles of the theory will be quickiy 
reviewed, but the study will focus on the problems that 
occur when dealing with nonlinear theory. Most of ths 
information has been taken from ~Ref. 6] as w2ll as from 
the course the author had at M.I.T. with K.J. Bathe in 1982. 


B. THE NEED FOR A NEW THEORY 


Considering a coordinate frame defined by (i,j,k). a 

body of volume Vin which point A(x1,%2,%,) is subjected to 
the displacements (Uj ,,U5 /U3), GOGEeCSpending €O a Str ain 
vector {fe} (figure £8.11). 
In the following sections, the superscripts 0 and ¢ will 
refer to the body at time 0 and t respectively, the 
Subscripts 0 and t will refer to the configuration at time 0 
and t respectively. This chapter, HOE Pens pusvose vor 
Sraplicity, will first deal with the static nonlinear anal- 
ysis of the material. 

In the linear theory of finita elements, one uses the 
well known Cauchy stress tensor 1; associated with the 


engineering strain tensor aj; “YS + = | : 
ox 1 
Then, the principle of virtual work 1s written: 
where *R represents the virtual work of externally applied 


Ferces and §6e€ is a virtual strain. 
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Figure B.1 Geometric Conventions. 


< 
| Bem * R (B.1) 


Equation (B.1) is then discretized over the body and becomes 
a set of integrations over each of th2 finite elements. In 
the case ofa large displacement, the volume of the body 
over which ths integration is performed might have signifi- 
cantly changed. Also notice that 2quation (B.1) is written 
in the original coordinate frame (diefined at *=)) and that 

Ty 2nd Gyn refer to the current configuration of the body. 
The Cauchy stresses at time t+#At cannot be obtained by 
adding an increment due to the straining of th2 material to 
the stresses at time t . The rigid body rotation of the 
Material has to be taken into azscount since the Cauchy 
stresses vary under rigid body motions. Therefore, we must 
perform the integration of equation (B.1) over the unknown 
current volume with respect *9 th2 original g2ometry that 
could be significantly different from the current one. 

The above discussion e2@mphasizes the need for a naw set 
of stress and strain tensors that would alleviate the afore- 
mentioned problems and enable the integration of the prin- 
Ciple of virtual work to be performed. 
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C. DEFINING HEW STRESS AND STRAIN TENSORS 


1. Green-Lagranage Strain Tensor 


The structure introduced aarlier is considered and 
two points A and B at t=0, or coordinates A (°x) P “B("X) 
are defined . At time t, the body has been deformed and A 
and B have moved to “A("x) and “B(*K) ($a guce Ba2). 


ee et ger ee eet ee ees ee ee 





Figure B.2 Displacements Conventions. 


A Taylor expansion is used to express the coordinates of B 


mea function of the coordinates cf A. 


Me = Sy + OE, Oxp XY (B.2) 
235 
eee With d *x.= oe = %, and d “x; = “x; > aa , gives : 
fee.) = (9%, )(4 x) (B 3) 
5 x 
meas; = «6 Gx) (a °x} (B.4) 
where [ ox] = ( 3x; ) ¢ (dete = (i “25 ¢ (a °x} = (a °x; ) 
ox} 


In other words, equation (B.4) expresses how a small fiber 


defined by the vector {d %} at ct=0 has rotated and 
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extended between time 0 and time t when it becomes (aé x} 
The matrix od is called the "deformation gradient" 

& ; 
The new length dS of the fiber will be : 


("ase = fax} ° (a°x} ( T refers to the transposed 
matrix) 
and therefore 
Ome 2 t 0 

meeps = (5x}t4 x}) (5% fa -x}) 

t 0 ec 0 

meaoye = {d x} (,C] (4 x} (B55) 
te i ed ; 
mee = [ort] {[oX] iss called the "Cauchy-Green deformation 
fomsor’, 


Notice that if the Cauchy-Green deformation tensor is iden- 
tity, equation (B.5) indicates that the length of any fiber 
will not vary. In other words, whenever rigid body motions 
are considered, the Cauchy Green deformation tensor is iden- 
meey sance the fibers do not stretch. 

The principles of the finits element method will now 
mepagercKkiyetecalled. Assume that the displacement (u,) of 
any point of the body can be written as a combination of the 


displacements of N selected points called "nodes" ;: 


(a) = acs) (B.6) 
k=1 


Where the h are interpolation ‘functions that depend only on 
th geometry of the body. In addition, the nodes are chosen 
so as to get a division into quadrilateral elements and it 
1s assumed that the displacement sf any point is only 2 
function of the displacements of th2 corner nodes of the 
element it belongs to. Then equation (B.6) becomes 


4 
(u; ) =) me OT (B.7) 
k=1 
Recalling that (u;) is Simply (x3) = (°x5) gJives: 
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L 


4 
(x) = (x) >. Been (B.8) 
K=1 


e ; : z 
The derorma*tion matrix te Xj] can be expressed very simply 


using the previous equation : 


($X] = @ R) = 4, -y dhe (a) (B.9) 
Ez k= dX5 


Define row the "Green Lagrange" strain tensor by : 
t t 
Mer = 172 ((,C) - (ry (B. 10) 


Where [I] is the identity matrix. 

From the previous derivations, it can be observed that the 
reen Lagrange strain tensor is ) for rigid body motions. 
The Green-Lagrange strain tensor refers to the body at «ime 
neem respect to the initial configuration. This is why it 


will be so useful in dealing with large deformations. 


Recall that the ultimate goal is to apply the prin- 
Ciple of virtual work to the structure under study. In 
particular, having defined a new strain tensor, the relation 
Giving the virtual Green-Lagrang? strain tensor corre- 
sponding to a virtual displacement (Su) must be known. Tn 
the case of a linear problem, the virtual engineering strain 
tensor would be: 
O25, = 1/2 ut Quy) 
In the case of large displacements, the Green-Lagrange 
strain tensor should be used. It is shown in [Ref. 6] that 


the virtual GSreen Lagrange strain tensor corresponding to 


Virtual exgineering strains 2S: 
$ 0%; 7 5X, ax, St Gn (B. 11) 
94 9x; 
or : 
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S25 = 8x, ox; 608i; (B. 12) 
re 


Having defined a strain tensor which is invariant under 
rigid body motions, a stress tensor corresponding to the 
Green-Lagrange strain tensor needs to be defined. 


2. Stress Measures 


Starting with the Cauchy stress tensor . she 
Piola-Kirchoff stress tensor is defined : 


c ( 
O° i; “= oe arn X5n (Bs43)) 


Where =, are the densities of the material at time 0 and 
t respectively and ( kag = an 1s the inverse of the 
deformation tensor defined previously. 

Equation B.13 can be easily rewritten in equivalent forn: 


Rt =< 5 xin Si Xj n (B. 14) 


It can also be shown by using the principle of mass conser- 


vation that : 


toe = to detf 5X] (B. 15) 


D. PRINCIPLE OF VIRTUAL WORK 


The principle of virtual work of equation (B.1) is 
rewritten using the strain and stress tensors defined previ- 
ously in equations (B.12) and (B.14) ; 


ie ey c c 
ae = 0X oS OMk © Mim EXP GH ei; dV (B. 16) 
G 


OT: 
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E eae fe ee 
R = 7 054,S0cif¥ (B.17) 
V 
and using equation (B.15) 


= 6Si;60ei AV (Be ic} 
Oy 


Thus, the principle of virtual work has been simply 
expressed in terms of a new set of stress and strain 


tensors, integrated over the origin2l volume ofa 


Eo. THE INCRESENTAL CONTINGUM MECHANICS EQUATIONS 


In this section, the principle of virtual work will be 
applied tothe structure and the incremental formulation 
using the Piola-Kirchoff and Green-Lagrange tensors wiil be 
developed. Non linear terms will arise from the rather 
complicated definiticn of the strain tansor, but it must be 
pointed out that the new formulation orovide the means for 
the modelling of large deformations. 


ct 
i" 


Assume that the configuration of the pody at time a 


known, the configuration at time t+At must be determined. 
Writing the principle of virtual work at time t+At gives : 


t+ At ome 
9 515 SCAT = R (B. 19} 
Orz 


tHR +: external work at time t+#At 
C, bt e *® ® e 
§ o€i3 Virtual increment 1n G.L. strain 


bet : stress state at time t+At 


Separating between the known terms that refer to the config- 
uration at tine + and the unknown terms that are the incre- 
ments of stress and strain betwean time t+ and time trAt 


gives : 
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t+ At £ tat c 
— 95 = 9 543 * OSij and — = ofijtofi; © (B. 20) 


Meeando&is are simply the incremants in stress and strain 
r2spectively. 

Te is derived in (Ref. 6] that the increment 0 
Green-Lagrange strain 0f&ij is made of a linear part oé¢@ijand a 
nonlinear onegl;. - The term lineac refers to the increment 
in displacement u; . 


OF = iz toh, cr © oy = 86 ( 044 tony ) (B. 21) 
Using equation (B.21) in equation (8.19) gives : 


[ o533 + Sid ( 8 (oeij to1ig)) ave = MOR (B. 22) 
O.7 


Again, separating between the known and unknown terms gives 


: 0Sif i dV | oS 8ohig = oar -| 9Sij 608i; (B. 23) 
Q7 Vv oy 
For small increments in displacenents, equation (8B. 23) 
written at time t+ indicates that 5 oeij= Soe 4; - This signi- 
fies that in equation (B.21) the aon Linear term is negli- 
gible 

Then the constitutive law of the material detailed in 


chapter II allows to relate stress2s and strains: 


054 = OGeshij ~ Cizs 0F% (B. 24) 


and B.23 becomes: 


[08 Cips sd V | 6 Sy SonhjdV= TR [sy i549 (B. 25) 
V ny 


Ov 
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Tie tight hand side terms of the previous equation ars 
known: 

(Pi; 15 the linear increment in virtaal strain involving only 
known ‘terms 

0Sij is known from the previous time step 

*ae is the work of external prescribed virtual forces. 

The left hand side of B.25 is unknown since o%jand Oy; 


involve the displacement from time t to time t+At . 


Feo FINITE ELEMENT DISCRETIZATION 


Equation (8.25) will be discretized over the structure, 
uSing the finite element approximation defined in section 
be . Let N be the number of naodes, the principle of 
Virtual work will be invoked N times, setting a unit 
displacement at each node in turn: 


6 u = 1 e Su, = 0 /KFj 


A system of equations whose unknowns are the nodal displace- 
ments is cbtained Let {Au} be the vector of unknown 
displacements, {PF} be the vector of nodal point forces 
equivalent to the internal stresses. 

Then equation (B.25) can b2 rewritt2n in matrix form : 


Doky )@u} + Ch,1fm3 = FS) - PR; (B. 26) 


em.) and C %&_) are known from the material characteristics 
at time + and correspond respectively to a linear and nonli- 
near contribution. It is therefor2 possible to solve equa- 
tion (B.25) for {Au}. Howaver, because of the assumption in 
equation (B.24), the exact solution might not’ be reached 
immediately. Furthermore, depending on the «ime step size, 
the solution process might even be unstabie! In any case, 
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an iteration for solving =quation (B.26) must be performed 
until che exact solution is reached . 
A widely used scheme is the "modified Newton iteration" 


defined by the following equation and boundary conditions : 


Mees eikeg) Our = YR - te i? (B. 27) 
and @o8 is wad (Ay) pede en ee ni 2 aliecsnditionsa«: 
Fy = cty and {Fy = {FP} 


{au}* is the vector of incremental nodal point displacements 
at iteration i. 
aA, ' 
{ R} is the vector of applied loads (cons*ant in the iterza- 
=O ) 
eh el ; 
{ F} is the vector of nodal point forces equivalent to the 
stresses at time t+#At, iteration i-1. 

At the first iteration, equation (B.27) reduces to equa- 
ie: , , i 
tion (B.26) giving an increment of displacement {Au} . Then 
a better approximaticn of oe, ° O04 1s obtained . 6545 is 
updated to the new state of stresses and becomes 9 5: : 
Equation (B.27) is then used to determine the new increment 
in displacement {Au} pe and) SO On Sun~1tt the -inerement in 


te Ab txt a 
{ 2} 


displacement is small enough, so that -eieetes) 17 


equation (B.27) 


G. INCLUSION OF DYNAMIC FORCES 


If the loads are applied rapidly, inertia forces need to 
be considered and a truly dynamic problem has to be solved. 
Using d'Alembart's principle, the element inertia forces are 
Simply included as part of the body forces. Let (u} be the 
vector of nedal accelerations and [M] be the mass matrix of 
the system. Then the principle of virtual work is written in 


the follcwing way 3; 
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te fe te AC s- ©6266 
Mee li au} + [ok jfAw} = RR -[N} Cau} - F (B. 28) 


Equation (8.28) represents a system of differantial squa- 
tions of second order. Gee is 990 Linear tezm (Kor wers 
negligible, the solution could be obtained by standard 
procedures for solving differential equations with constant 
coefficients. However, the procedur2s proposed for the solu- 
tion of general systems of differential equations can become 
very expensive if the order of th? matrices is large. 
Therefore, whenever the system is linear or nonlinear, scme 
effective methods for solving ‘the equations governing the 


equilibrium are required. 
1. Direct Integration Methods 


The essence of direct integration methods is based 
on two ideas. Estee 2S) almed. to satisfy B.28 only at 
Seeeain time intervals apart. Second, a variation of accel- 
eration velocities and displacements is assumed within each 
time step. The form of the assumption determines the accu- 
macy, stability and cost of each method. 

In the following, assume that the initial conditions 
(displacements, accelerations, YVerocitz@es) at tine 0, 
denoted (% a aa ) are known. In the solution, the tine 
Span under consideration, T , is subdivided into n equal 
time intervals At. Assuming that the sclution is known at 
time t , the methods of getting the solution at time t+At 
will be investigated. 


2. Central Difference Method 


In the Central Difference method, a finite differ- 


ence approximation will give the acceleration at time t ; 
- rAt t E 
‘y= 1 al ae ee th) (B. 29) 
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mereang the principle of virtual work at time <* and 


Sm@eetaituting lnto equation (B.30) gives then: 


cujeéuy + (5K Wf uy = ctr} (B. 30) 


t+ t te eo: tt 
Mees = (BY - Mok] - Amity Cay Baye uy (8.34) 


The previous equation gives the deformation at time ttAt 
from the characteristics of th2 system at time ct. 

When [M] is diagonal, which is frequently the case 
for mass matrices, tae sOmlc2 Og jae  tzme Et At does not 
mavolve any triangular factorization of the matrix [M4] , 
thus leading to more efficient computations. 

The shortcoming in the use of the central differance 
method lies in the time step restriction: for stability, thé 
time step siz? t must be smaller than a critical time step 
Been ich is equal to fT / , where T is the smallest period 
cf the finite element system. 


The central difference sch2me is fairly easy <*o 
implement for the integration of a system of non linear 
differential equations. However, because of thea limitations 
of the cime step , it might not be suitable for cases when 


loads are varying at a slow pace. 


3. Implicit Integration Schemes 





Since the interest of this study lies in central 
difference schemes, this section will be limited to a short 
description of the fundamentals of implicit time integration 


schemes. 
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Implicit time integration schemes use the princirle 
Meevaitrtual work written at time t+ At and not at time t as 
for the central difference method : 


Ab tit tt trAt | 
CH) {uy + (K] {vu} = ( &} (3. 32) 


Again, using a finite difference approximation of Sr and 
replacing into equation (8.32) enables to solve for { i} 
Since the formulation involves the rigidity matrix {6 Kj and 
the external work eR which are both unknown, the system 
has to be solved in a Similiar way to the Moditied Newton 


iteration that was detailed previously. 


Implicit time integration sthames are stable regard- 
less of the size of the time step used. However, if the time 
step size is too large, significant 2rrors can be accunu- 
lated a* each time step, leading to unrealistic results. 

The reader will find more details on the various 
miplicit method in [ Ref. 6]. Yo Lt ce Gampoint ed Cut that 
implicit methods are more tedious t9 implement. On the other 
hand, a larg2r time step can be used in the solution proce- 
dure, which can be of extreme importance when tudying 


phenomena over a Significant period of time. 
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APPENDIX C 
HOW TO USE THE TRANSLATOR EPSA-PATRAN-G 


This appendix is intended to 2xplain in detail the use 
of EPSA/PATRAN-G post-processing facilities. It is divided 
in three sections that will deal succ2ssively with 1) the 
displaying of the original model ; 2)czhe deformed geometry; 


3) the cortour plots of elament and nodal points quantities. 


A. DISPLAYING THE ORIGINAL MODEL 


When making an initial EPSA analysis ona particular 
structure, the geometry of the moiel has to be input inte 
PATRAN-G. As explained in chapter [II, all the geometrical 
information is contained ina file FORO19.DAT that is 
created each time an EPSA run is nade. rie Lhouc, Of she 
Original geometry must be made via the neutral input mode of 
EPSA. The procedure, starting fron the "logon" to PATRAN-G 
Meeche £cllowing : 

- Select the SO option 

- Select the new data file option (option 1) 
- Select the neutral system (option 4) 

- Select the input mode (option 4) 

- Input the neutal file name : FORO19.DAT 


The original geometry will then be displayed on the screen 
It is often found convenient to have a perpective view 

of the model under study. In that cas3, the user should : 
- Issue the VIEW command 
- Select the rotation about the absolute axes (option 1) 
- Input an angle of rotation 

(23,-34,0 will give a very nice view but any angle can be 
input) 
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- Issue a PLOT command to have tha model displayed in the 


new axes. 


An example of the procedure is provided on table VI. 


: 
| 
a 


TABLE VI 
Finite Element Model Input Procedure 


MODE? 1.GEOMETRY MODEL e.ANALYSIS MODEL 3.DISPLAY 4.NEUTRAL SYS. S.END 
4 


NEUTRAL FILE? 1.CREATE OUTPUT @.INPUT MODEL 3.POST-PROCESSING 4.END 
de 
INPUT NEUTRAL FILE NAME 
2FOR@19.DAT 
DO YOU WISH TO OFFSET ANY NEUTRAL INPUT IDS? (Y/N) 
aH 
EPM e0@ STEPS ,NO STIFFENERS W=3@. 
SHALL WE PROCEED WITH THE READING OF THIS FILE? (Y/N) 
y 


Se eee se opin ee ee OM, coo) ee) Se es i ge ee SS eee ee ee ee ss Se 


The PLOT command can be issued anytime to display the orig- 
inal geometry on the screen. 

When studying ccmplex models, one does not want the 
element and node numbers *t9 be printed along with the geon- 
etry of the structure. The command SET, LABS3, OFF followed 
by a PLOT will display the original geometry withou+ any 
labels printed. 

When a model has been input into PATRAN-G, a data file 
PATRAN.DAT is created on the user's directory. When 
connecting with PATRAN-G at a later time, the user can 
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meeeerce the option "“Jast data file " (option3) to have the 
original geometry displayed on the screen without having to 


input the model again. 


Be DEFORMED GEOMETRY 


On «the second card of the PATRAN-G input deck, the user 
specifies the number of PATRAN-G displacement data file and 
the number of element results data fil2. Two non-zero inte- 
gers in free format must be placed at the end of the second 
card (secticn II in the user's manual) requesting the number 
of displacements and of elements files respectively. 

Assuming that the user has mad2 a 200 steps run with 10 
output requests for PATRAN-G displacement files, ten (10) 
files FOROI8.DAT will then be created at equal time inter- 
vals. The deformed geometry corresponding to time st2p 100 
will therefore be contained in FORO18.DAT;5. 

To display the deformed geometry corresponding +o time 
step 100, the user should issue the following commands: 


RUN,DEF : requests deformed geometry option 
Input the name of the displacements file :FORO18.DAT 5 
Semect the PLOT option (option 3) 


- Select the undefcrmed gecmetry (2) followed by the 
deformed geometry (3). An example of the procedure is 
provided on table VII. The undeformed geometry superposed 
With the deformed gecmetry will then be displayed on «he 
screen. 


C. POST-PROCESSING OF ANALYIS RESOLTS. 


Element-related quantities 1ik2 von Mises stresses are 
contained in FORO16.DAT files, noial point quantities are 
stored in FOROQ18.DAT files. As described in chapter III, 
each column of those fil@s contains a specific quantity 
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TABLE VIL 
Deformed Geometry Procedure 


ae oe eee tee a 


> 
MODE? 1.GEOMETRY MODEL 2.ANALYSIS MODEL 3.DISPLAY 4.NEUTRAL SYS. S.END 
>RUN,CON,COL, 4 
CURRENT FILE FOR NODAL RESULTS IS FOROQ18.DAT 
NEUTRAL RESULTS FILE? 1.NEWU FILE e.CURRENT FILE 


1 
INPUT THE RESULTS FILE NAME: 
FORO18.DAT;S 

DATA WIDTH = S 


Pace FITLE =EPM 800 STEPS 
7.5000129E-04 


DATA VALUES RANGE FROM -0.5e5E+00 TO @.134E+00 
ASSIGNMENT? 1.AUTO @2@.MANUAL 3.SEMI-AUTO 4.USE CURRENT LEYELS 5.END 
»4 

PREVIOUS CONTOUR LEVELS USED. 

MODE? 1.GEOMETRY MODEL e@.ANALYSIS MODEL 3.DISPLAY 4.NEUTRAL SYS. 5.END 
»RUN,HI,C 


i  ecciaeeini ——— cy I ag a I ct ll ee I mess a es igs le, 


eg, cee cement ome Oe ete ee ey ee et ee ey SE GES go SE ee ee 


| 





Ge. 2. column 31 of FORO16.DAT contains the von Mises 
stresses). The reader will refer to chapter III for the 
detailed organization of those files. 

Again, assume a 200 time steps analysis, with 10 output 
requests for element results files. The user might want to 
display the contour plots of von Mises stresses 2t time step 


100 . In this cas# the following commands should be input : 


- RON, CON, COL, 31: tells PATRAN-G to look at column 31 
that contains the von Mises stresses. 

- Input the file name FORO16.DAT 5 

- PATRAN-G will then ask for a color assignment (automatic, 
manual, semi-automatic, current levels used) that the user 
will select according to his needs. 

The contour plots are then ready to be displayed : 

SenUN, Hi, FR would display the color-filled contour plots 

- RUN, HI, CON would display the contour plots (color lines) 
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The von Mises stresses are the nost useful elemen= quan- 
tities to be displayed, but other elemant-relatesd quantities 
detailed in chapter ITI like the x and y stresses in local 
or global ccordinates could be displayed as well by looking 
at their corresponding column in the element results data 
file. 


Dealing with nodal point quantities, the displacement 
normal and the velocity normal to the shell are very nean- 
ingful quantities in an analysis. They are respectively 
stored in column 4 and 5 of the displacement results files 
FORO18.DAT . 

A contour plot of the normal displacement at time step 100 


would then be obtained via the following commands: 


RUN, CON, COL,4 (look at column 4) 
FORO18.DAT 5 (name of the £1182) 

- Color assignment chcsen 

RUN, HI, C or RUN, HL, f&R 


' 


Notice Woeca the ftiwid-structure interaction is ON, the 
normal displacement contained in csolumn 4 corresponds to 
pure deformations, the rigid body contribution having been 


taken out. 


All the element results processing is implemented in the 
routine NEUSTRE, alithe nodal points processing is imple- 
mented in routine NEUDISP. It should be pointed out that any 
modification to the capabilities of the translator (1.2. 
being able t¢)9 display other types of quantities) can be 
made by modifying these routines only. 
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APPENDIX D 


PP 
LISTINGS 

This appendix contains the varisus files that constituts 
the translator module. The submodul2 that displays the orig- 
inal geometry is listed on the first four pages. It is 
imbedded in file FRANK.FOR. The subnaodule thax takes care of 
the post-processing facilities is imbedded in file DISP.FOR 
and is listed in the remaining pages. It has been mentioned 
previously that EPSA had been slightly modifisd to accome- 
date color graphics capabilities. Thewoniy Interaction “Of 
EPSA with the translator occurs vl3a subroutine calls. All 
mee “calls" occur in COMPUTE (f5r post-processing) and 
REPORT (for the original geometry). A labelled COMMON 
called FRANK has been created and is défined in the routine 
AAA as reguir2d by EPSA. The requests for PATRAN-G outputs 
are echoed in the EPSA output £113, all the nodifications 
for that purpose having been made in routine READIN. The 
user must be aware that the size of blank COMMON array A has 
been increased to store the deflections in x, y, 2z direc- 
tions instead of only the z direction previously. 
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AAA A 


n 


10 


11 


ie 


SubDroUtIne orelim 
dimension A(1) 


common 


1a(1) 


eauivalence (Cia(l),al(l)) 


common/ssize/ 


noauad, 
nsstyor 


norot 
nntot, | 


common /cPara/ nsteos,s 


common /stab / 


ioe C1). j 


iba(1l)snajr,-neltot,nibd,nioadsnbrect, 
isheet, 
AN je 


Se nsotse 
gcdso, liaut, Ibcale 

noegr, nends nsheet 
Ssize,jsoar,jivelo,jstrersjxmas,jielm, jomat, 


nstrotsr nvotse nhpts, 


jllods,jloco-vjoret,jlhiss,jistrn,jforcrejxelocejnair,janis, jnadeg, 


jlside 


common/frank/nfntotsliu 
COMMON/TITRE/NTITLE(80) 


CHARACTER#9 
CHARACTER#S 
CHARACTER*#8 
CHARACTER*9 


this routine comoutes new number of nodal 


for any 


BUFF 
CUTE 
TIM 
VER 


given sheet 


KJSJNGBEG=JNNI 
nfntot=nntotelTACJINNI) fA CINGBEG=1) -2*(KJ=2) 


Tltypese 
Ylkesl - 
}livsd 
T}iasd 
lini=9 
Tin2=0 
Tin3s0 
Tin&s0 
}inSs0 


weite(llu,19) 


Titype, Id) 


format lie, 318) 


write(liu-11) 


Cenc Cl 


format(80 Al) 


points 


idellivelilkes,lIinislined,tin3d-lIn4, Vind 


),121,80) 


taking care of second oacket 


litype=e6 
likc=l 
Jinlenfntot 
lind=neltot 


write(llu,10) 
dJate(BuUFF) 


call 


Tityoe,1) 


write(llusle) BUFF 
format(A,3xX%,'17:122909',5%X%,'1.4") 


return 
end 


subroutine shconn 
dimension A(1) 


common 


1a(1) 


idsllivellke,linizlines,lTin3d, Ing, ins 


equivalence (Cial(l)-,a(l1)) 


common/ssize/ 


noaquad,s 
Nsstvor 


1Sheet, 
NNie 


iball),na 
nornoet 
nntot, | 


jeneltotsnivdas,niloadsnorect, 
Se NSOtser nstrotss, nvetse nhpets, 
gaso, li90g, locale 


oa) 





aaqana 


aAanNgn 


a 


common /codara/s nstensr nbeyjs nends, nsheet 

common /stab / ibtll),jssi2ze,jsoar,jvelors,j;strers,jxmase,jielm, jomat, 
i jlicdejlodo-,joret,jilhisrejsternsiforer,jxloce. jnqisjnnise jnaneg, 

2 j lside 


common/frank/nfintotsl lu 


this routine will write the element corner nodes of each 
element on neutral file 


lltype=e 
llkese 
linod=4 
llive4 
LLN1=0 
LLN2=0 
LLN3=0 
LLN4=0 
LLNS=0 
LLCONF =0 
LLPIO=0 
LLCEI=0 
THET1=0. 
THET?2=0. 
THET3=0,. 


INITIALIZATION DONE 


noreyv=0 
llelsjnniejnai 


do 200 k=l, llel 
llerow=ITACjnaitk=1) 
nmunder of elts in each row 
do 100 j;=1,llrow 
lel=IA( jnapeg+Kel)+j-1 
nlell=j+noreyv 
nlel2=jtitnorev 
nlel3sjtltnorevtllrowet! 
nlel4=jtnorevetl lrowtl 
ready to disolay oacket 


JIid=LlEL 
weire(llu,80) lityoe,llidszllivellke,-LLNI,LLN2,LLN3,LLNG,LLNS 
80 formatl(i2,8i8) 
weueecl)lts,s!) |limod,-LLCONF,LLPIO,LLCEI,THETI,THET2, THETS 
81 format(i18,318,3e16.9) 
weite(llus,82) nlellsniel2,niel3rsniel4 
82 formart(10i18) 
100 continue 


norev=norevetl lrowel 
200 continue 

return 

end 


Subroutine sheet f 
dimension A(1) 

common jal) 

equivalence (Cia(l)-a(l)) 
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an 


ANN 


o 


commnon/ssize/ togl(l),-najyrneltots,nidot-nloadsnbdrect, 

l Noaquad, isheet, Noratsse NSOCTS, Nstrotssr Nvotse NNOtS, 

EeeemsSst vor, ANMj+s Antots I39dsor liqud, tbcate 

common /coara/ nst#osr nbers, nends, nsneet 

common /stan / itt(l)-jssizerjsnar, jvelo-,jstresjxnass,jielm,sjomat, 
1 jllinods-jlodo,joret,jlhisr-jstran-jforec,jxlocre jnairsjnnisjnaveg, 

2 jlsice 


common/frank/nfntotellu 
character*l gtyoe 


renumbers the nodesr chanae coordinatess disolay oacket 1 


lityoe=1 : 
Ttida=O 
Ttiv=s0 
llke=e 
Tintsd 
ltn2=0 
}1n3=0 
T1n4ds0 
11n5=0 
ncoun=0 
ltdof=6 
lliecf=l 
atyoe='G' 
lliconf=0 
Tleidg=9 
LSPC1=0 
LSPC2=0 
LSPC3=9 
LSPC4=0 
LSFCS=0 
LSPC6=0 


INITIALIZATION OONE 


kj=jnabege-jnni 
krowskj-2 
do 200 j=l,krow 
loom on new nbe of row 


neasunszncountlIACjnnitjel) 
TloteIACjnnitjeli-e2 
do 100 l=t,ltor 
Tlid=llidtt 
xlloczACncoun*2tjxloctre*!) 
ylloc=A(necounxetjxloct2xt +l) 
skio first node of row 


zitoc=0. 

if (ACjvelow2d).ne.9.) then 
reouro=l1./A(jvelosze) 
theta=xlloc/rcouro 
xlloc=rcourbresin(theta) 
ziloc=rcourbrcos(theta)=rcourbd 
endif 


if (a(jvelo-l).ne.0.) then 


reouro=l.e/aljvelonm-l) 
theta=viloc/rcourys 


57 





70 


a 


80 


81 


82 
100 
200 


80 


ylloe=rcouror*sin(theta) 
zlloc=rcourorcos(theta)orcours 
endi f 


if€CalCjvelowrl).ne.0.).and.lCaljvelow-2).ne.0.)) then 
weite(llu,70) 

format('error,two curvatures are non zero’) 
endi f 


ready to disolay voacket 


weite(]1u,80) Tltyoe,llidszllivellke,linistineds-1in3-,}1n4,11905 
format(i2,8i3) 
write(llu,-81) xllocrs,yvylloc-ziltoe 
format(3e16.9) 
write(llu,82) llicfrsgqtyoer,lidofsllconfrellcids-LSPC1,LSPCe, 
moPecs LSPC4,LSPCS,LSPC6 
format(I1,a1,318+2x.6i1) 
continue 
continue 
return 
end 


subroutine endaneu 
common/frank/nfntot, I iu 
litype=99 
Tlid=90 
Thived 
like=l 
weite(llu,80) iltyoe,llidsllivelike 
format(i2, 318) 
return 
end 
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SUBROUTINE NEVOISP(DEFL,VELO) REPORT 


OLMENSION A(1) BLANC 2 

COMMON IA(C(1) BLANC 3 

EQUIVALENCE (IAC 1), ACI) ) BLANC 4 

BEIMENSION DEFL(1)+VELOC1) 

COMMON /CPARA/ NSTEPS , NBEG , NENO , NSHEET - N2BO , N3BD1 , CPARA 2 
1 N3B802, INTRVL, OELT, NHTOT, NJOIN, NRELAX, ALPHA CPARA 3 
§$ LEN SBX . MSPARA 3 

COMMON / SSIZE / I8G6(1)-/NOJ,NELTOT,-N1BO-NLO40,NBRECT, SS UZE > 2 
1 NBGUAD, ISHEET, NPRPTS, NSPTS, NSTRPTS- NVPTS, NHPTS, SSIZE 3 
2 NSSTYP, NNJ, NNTOT, LGOSP, LIQUO, LBCALC Samee, 4 

c REPORT 8 

COMMON /STAB / IBT(C1)-JSSIZE,JSPAR,JVELO-JSTRE,JXMAS,JIELM,-JB™AT, STAB 2 
1 JL180,JLOOP, JPRET,JLHIS,JSTRN, JFORC,JXLOC,JNGI,JNNI,-JNQBEG, STAB 5 
$ JLSIDE-JIELMCL-JSTIF,JOEFL,JFORLG, STAB 4 
2 JIFPAR, JFLPAR, JXCOORO, JYCOORO,JOELTAX, JOELTAY,JVMA,JSEFX, STAB 5 
5 JFLUFR, JPRINC, JVELRAD, JGENFR,JPRES,JCSEP STAB 6 


COMMON/CIO/ NIN,NOUT,NTHIST,NCORT,MCORT,NTIPLOT,NVMA4,ITITLE (20) 
COMMON/FRANK/NFNTOT,LLU-LLN,LLS,NOPLTS,NSPLITS 
COMMON/TITRE/NTITLE (80) 

GomMON/COELT/ISTEP, TIME 


aan 


INTEGER NSUBT1(30) -NSUBT2(80) 


KJSJNGBEG-JNNI 
NFNTOTSNNTOT*IACINNI) LAC JNGBEG=1) -2%(KJ=-2) 
LLU=19 
LLN=18 
LLM=17 
DEFMAX=0. 
MAXNOD=SNFNTOT 
NWIOTHSS 
NOMAX=0) 
LLIO=0 
NCOUN=0 
KROWSX Jee 
00 10 J=1,80 
NSUBT1(J)=0 
10 NSUBT2(J)=0 


Peure tHE TIME STEP ON FIRST SUBTITLE: 


CLOSE (UNI T=LL‘™) 

OPENCUNTIT=LLM,STATUS='NEW® ) 

WRITECLLM,*) TIME 

CLOSE CUNIT=LLM) 

REAOQ(LLY,99) (NSUBTI(J)-J=1,80) 

CLOSE CUNIT=LL™) 

OPENCUNTT=LLM,STATUS='OLO') 
OPENCUNIT=LLN,FORMsS'UNFORMATTED',STATUS='NEW! ) 


KROW IS THE NEW NUMB. OF ROWS ,KJ IS THE OLO 
NCOUN COUNTS NODAL PTS IN EPSA, LLID COUNTS FOR ACTUAL MOOEL 
NUMR COUNTS NOOAL PTS IN EPSA 


FIRST LOOP TO GET MAX, OEFORMATION AND NOOE NUMBER 


aAaNnNnNANANnNG 


00 200 J=1,KROW 


a3 





AANnANnNnANHn 


AAO 


AA 


aan 


Aan 


100 
200 


90 
95 


91 
99 


NCOUN=NCOUN4+ TAC JINNI+J=1) 
LLPTSIACJNNI+J) =e 


LLPT IS NUMS. OF POINTS (FOR PATRAN) IN EACH ROW 

we TAKE THE RIGID 800Y YOTION OUT BY SUBTRACTING THE V AND w 
DISPLACEMENTS AT ENDSLOCK AT EACH NODAL POIPIT 

NCOUNt2 IS THE POINT NB. CEPSA) FOR ENOSLOCK 


RBYZDEFL(3*(NCDUN+t+2) 1) 

RBZ=DEFL(3*«(NCOUN+2) ) 

[PE ScCET GUO LE. G) — THEN 

REY¥=0. 

R6Z=0. 

END GF 

DOMmTO 0s L=1,L LPT 

LLID=LLID+! 

NUMR=NCOUN+L?# 1 
OXSDEFL (C3 *NUMR=-2) 
DYSDEFL(3*NUMR=-1)-RBY 
DZ=DEFL (3 *NUYR) RBZ 


IS SHELL CURVED,CHANGE COORDINATES 


Meeeea tJ VELD=1).NE.0.)-0R.CACIVELO<-2) .NE.9.)) THEN 
CALL CHCODROCNUMR,0xX,D0Y,0Z) 
ENOIF 
DD=AMAXICABS( OX), ABS(OY)-ABS(0Z)) 
PR CASS (DEP MAX) .LIT.D00) THEN 
NOMAX=LLIO 
DEFMAX=DD 
IF CABS CAMINI(0X,DY,0Z)).EG.00) THEN 
DEF WMAX=-DEFMAX 
ENOIF : 


WE CHECKED IF DEFMAX WAS NEG, 


ENOIF 
CONTINUE 
CONTINUE 
OK FOR FIRST CARD 
WRITECLLN) C(NTITLECTI),-121,80),NFNIOT,MAXNDD,DEFMAX,NDMAX, 


1 NWIOTH 


wreite(LL’,90) TITLE-NFNTOT,MAXNODO,DEFMAX,NOMAX,NWIODTH 
FORMAT(A,218,£16.9,218) 

FORMAT(2041,218,E16.9,218) 

WRITECLLN) (NSUBTICI),-I=1,80) 

WRITECLLN) (NSUBT2CT), 1=1,80) 

WRITECLLM,91) 

FORMATC'PINE') 

FORMAT(80A1) 

LLID=0 

NCOUN=0 


SECOND LOOP TD PICK UP DEFLECTION AT EACH NODE (COF PATRAN MOQDEL) 
DD 400 J=1,KROW 

NCOUN=NCOUN+I AC UNNI +J=1) 

LLPT=IACJNNI+J) <2 


TAKE OUT RIGID B80DY MOTION, RSBX,RBZ OISPLACEMENTS AT END BLOCKS 
ASSUMED 1D REPRESENT RIGID BODY “DITIONS 
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QRBYSOEFL (3% (NCOUN+2) 1) 
RBZ=OEPL( 32% (NCOUN+2)) 
Pe CE TOU0.E9.9) THEN 
RBY=0. 

RAZ=0. 

ENOIF 


C IF NO FLUIO OPTION DOO NOT SUBTACT RIGIO SODY CONTRIBUTION 


300 
400 


00 300 L=1,LLPT 
LLIO=LLIO+¢l 
NUMRENCOUN+L +1 
OX=DEFL(3*NUMR=2) 
DY=DEFL(3*NUMR=1) -RBY 
DZ=O0EFL (3 *NUMR) -RBZ 
S7EOC =0Z 

VELZSVELO( 3*NUMR ) 


Poca IVElO=1)) NE 00.) OR (AC IVELO=2) .NE.0.)) THEN 


CALL CHCOORO(NUMR,OX,0Y,02Z) 
ENOIF 
ARITEC(LLN) LLIO,0X,0Y,0Z,0ZL0C,VELZ 
ARTTEC(EE4, 72) LLIO,0X,0Y,02Z,0ZLOC,VELZ 
CONTINUE 
CONTINUE 


PeerOoe FILE OPENED ON UNIT LLN 


ann 


anNaAannn 
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19 


l 
S 


1 
2 


l 
§ 
eC 
3 


Gale CLOSECLLNV) 
FORMAT(18,5E16.9) 
RETURN 

ENO 


SUBROUTINE CHCOORD(N,X,Y,2Z) 

OIMENSION A(1) i 

COMMON IA(1) 

EQUIVALENCE (IAC 1), AC1) ) 

COMMON /CPARA/ NSTEPS , NBEG , NENO , NSHEET , 
N3802, INTRVL, OELT, NHTOT, NJOIN, NRELAX, 

weN SBXx 


N280 , N3B01 , 
ALPHA 


COMMON / SSIZE / IBG(1),NQJ,NELTOT,NIBO,NLOAO,NBRECT, 
NBQUAD, ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, 


NSSTYP, NNJ, NNTOT, LGOSP, LIQUD, LBCALC 


COMMON /STAB / IBT(1)-JSSIZE,-JSPAR,JVELO,-JSTRE,JXMAS,JIELM, JBMAT, 


JL180,JL9O0P,JPRET,JLHIS,JSTRN, JFORC,JXLOC,JNGI,JNNI,JNGBEG, 


JLSIOE,JIELMCL,JSTIF,JOEFL,JFORLG, 
JIFPAR, JFLPAR, JXCOORO, JYCOORD, JOELT4X,JOELT 
JFLUFR, JPRINC, JVELR2AD, JGENFR, JPRES, JCSEP 
COMMON/CIO/ NIN-NOUT,NTHIST,NCORT,MCORT,NTIPLOT 
COMMON/FRANK/NFEFNTOT,LLU,LLN-LLS,NOPLTS,NSPLTS 
OIMENSION XMAT(3,3) 
00 10 T=1,-3 
00 10 J=1t,3 
XMAT(CI,J)=0. 
CONTINUE 


AY, JVMA,JSEFX, 


oNVMA,ITITLE (20) 


lilo OURTNESCHANGES THE O1SPE. OF NOOAL PT. NC(FOR EPSA) 


IN A GLOBAL RECTANGULAR SYSTEM 


DPESGURVSTURE IN ¥ OIR. IS NON ZERO ,CALCULATE 


THE ROTATION MATRIX AT EACH POINT 


10 1 


REPORT 
BLANC 
BLANC 
BLANC 
CPARA 
CPARA 
MSPARA 
a9 lee 
So P2e 
So 1ZE 
REPORT 
STAB 
STAB 
STAB 
ST4B 
STAB 
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AAMN 


AaANaNM 


Pec ACIVELO-1).NE.0.) THEN 
RCQURY=S1./A(JVELO@1) 
THETALZACITXLOC +2*N-1)/2COURY 
XMAT(L,1LI=A1. 

XMAT( 2-2) =COS(THETA1 ) 

XMAT( 2-3) =SINC(THETAL) 

XMAT( 3,3) 2XMAT (2,2) 

XMAT( 3,2) =~1 .eXMAT( 2,3) 

CALL PROD(CXMAT,X,Y,2Z) 
ENOIF 


PeeCcURVATURE IN X DIRECTION IS NON ZERO: 


PeeC SC JVELO]2).NE.0.) THEN 
00 20 [=1,3 
DO 20 J=1,3 
XMAT(I,J)=0. 

CONTINUE 
RCOURX=1./A(JVELO=-2) 
THETAZ=A( JXLOC +2*N-2)/RCOURX 
XMAT(2,2) 51. 
XMAT(1,1)=COS(THETA2) 
XMAT(3,3)=XMAT(1,1) 

XMAT( 3,1) ="1.*SIN( THETA) 
XMAT( 1,3) =SSIN(THETA2) 
CALL PROD(XMAT,X,Y,2Z) 

ENDIF 

RETURN 

END 


SUBROUTINE PROD(XMAT,X,Y,Z) 
DIMENSION XMAT(3,3) 


THIS ROUTINE DOES THE MATRIX PRODUCT XMAT*(X,Y,2Z) 


X1=x 

Y1=yY 

Z1=2 

X=XMAT(1,1) *#X14XMAT(1L,2) *Y1+XMAT(1,3) 21 
YEXMAT( 2,1) eX1L+XMAT( 2,2) 241 4+XMAT( 2,3) #21 
Z=XMAT( 3,1) *#X1+XMAT(3,2) "71 +XMAT( 3,3) 21 
RETURN 

END 


SUBROUTINE PLOTOISP(DOEFL,VELO) 
BYMENSION DEFL(1),VELO(1) 
COMMON/COELT/ISTEP, TIME 


COMMON /CPARA/ NSTEPS , NBEG , NEND , NSHEET , N2BD , N3B01 , CPARA 2 
1 N3B02, INTRVL, DELT, NHTOT, NJOIN, NRELAX, ALPHA CPARA 43 
$3 LEN SBXxX MSPARA 3 


COMMON/FRANK/NFNTOT,LLU-LLN,LLS,NOPLTS,/NSPLTS 


CHECKS IF OUTPUT FOR PATRAN IS REQUESTED BY USER IF YES CALL 
NEUDISP 


fPeCiIiMe.EG.OEET) KDISP=1 
TIMEIT=FLOAT(NSTEPS) /FLOAT(NDPLTS) 
TIMEC=KOISP*#TIMET 
ITIMESJNINTCTIMEC) 
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anrany 


PeMeEC=TIME/CITIME*DELT ) 

ENOTIM=NSTEPS*DELT 
TESTSABS(TIMECeJNINTCTIMEC)) 
fieeetest.LE.0.00001).0R.(TIYE .EQ.ENOTIM)) THEN 
Paeey NEYDISP(OEFL, VELO) 

KOISP=KDISP+1 

ENOIF 

RETURN 

ENO 


SUBROUTINE PLOTSTRE(CSTRE) 

DIMENSION STRE(1) 

COMMON/COELT/ISTEP,TIME 

COMMON /CPARA/ NSTEPS , NBEG , NEND , NSHEET , N2BO , N38D1 , 
1 N38N2, INTRVL, OELT, NHTOT, NJOIN, NRELAX, ALPHA 

Ss FEN SBX 

COMMON/FRANK/NFNTOT,LLU,LLN,-LLS,NOPLTS,NSPLTS 


EmecKo [F OUTPUT FOR PATRAN IS REQUESTEO BY USER IF YES CALL 
NEUDISP 


PEeS(TIME.EG.OELT) KSTRE=S1 
TIMET=FLOAT(NSTEPS)/FLOAT(NSPLTS) 
TIMEC=KSTRE*«TIMET 
ITIMESJNINT(TIMEC) 
femee=TIME/(ITIME*OELT) 
ENOTIM=NSTEPS*DELT 
TESTSABSC(TIMECeJNINTCTIMEC) } 
ieeereo? .LE.0.00001) .OR.(TIME.EQ.ENOTIM)) THEN 
CALL NEUSTRECSTRE) 

KSTRE=KSTRE+1 

ENOIF 

RETURN 

END 


SUBROUTINE NEUSTRE(STRE) 

OIMENSION A(1) 

COMMON TA(1) 

EQUIVALENCE (ITA(1), A(1) ) 

OIMENSION OISP(10),STRA(10),0STRE(10),08S(5),VEL(2) 
OIMENSION STRE(1) 

COMMON /CPARA/ NSTEPS , NBEG , NEND , NSHEET , N28BO0 , N3B01 , 
1 N38D2, INTRVL, DELT, NHTOT, NJOIN, NRELAX, ALPHA 

Sy LEN SBX 

COMMON / SSIZE / I18G(1),NQJ,NELTOT,N1IB8O0,NLOAOQ,NBRECT, 

1 NBQUAO, ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, 
2 NSSTYP, NNJ, NNTOT, LGOSP, LIGUO, LBCALC 


COMMON /STAB / I8BT(1),JSSIZE,JSPAR, JVELO,JSTRE, JXMAS, JIELM, JBMAT, 
1 JLIBD, JLOOP,JPRET,JLHIS,JSTRN, JFORC, JXLOC, JINGI,JNNI,JNQBEG, 

$ JEOVOE,JIELMCE, Jot le, JOEFL, JFORLG, 
e JIFPAR, JFLPAR, JXCOORO, JYCOORD, JOELTAX, JOELTAY, JVMA,JSEFX, 

3 JFLUFR, JPRINC, JVELRAO, JGENFR,JPRES, JCSEP 

COMMON/CIO/ NIN, NOUT,NTHIST,NCORT,MCORT,NTPLOT,NVMA,ITITLE(20) 
COMMON/FRANK/NFNTOT,LLU,-LLN-LLS,NOPLTS,NSPLTS 
COMMON/TITRE/NTITLE (80) 

COMMON/COELT/ISTEP, TIME 
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CPARA 
CPARA 
MSPARA 


REPORT 
BLANC 
BLANC 
BLANC 


CPARA 
CPARA 
MSPARA 
So LZe 
SSIZE 
SSIZE 
REPORT 
STAB 
STAB 
STAS8 
STAB 
STAB 


e 


5 
3 


&wWwivt 
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COMMON/SPARA/HEG(1)-E-GNU,RHO,SGYLO,THICK,CURV(2) 
INTEGER NSUS8T1(080),NSUBT2(80) 
NAIOTH=31 
LLS=16 
ez=t5 
oo 10 J=1,80 
NSUBT2(J)=0 


G7TAKE CARE OF THE 2NO FITLE » INOICAT. TIME STEP 


C 


AAANAN 


oO AANAAAANAANAA 


801 


WE A 


O0S(4 


GEOSECUNTT=LEZ) 

OPENCUNIT=LLZ,STATUS='NEW? ) 

WRITECLLZ,*) TIME 

CEOSE CUNTT SLLZ) 

REAO(LLZ,801) (NSUBT1(1T),1=1,80) 

CLOSE CUNIT=LLZ) 

OPENCUNIT=LLZ,STATUS='OLD") 

FORMAT(80A1) 
OPEN(CUNIT=LLS,FORM='UNFORMATTEO',STATUS='NEW' ) 


RE GOING TO OISPLAY THE FIRST THREE RECOROS (TITLES) 


WRITECLLS) (NTITLE(I),1T=1,80),NWIOTH 
WRITECLLS) (NSUBT1(1I),1=1,80) 
WRITECLLS) (NSUBT2(1),f=1,80) 


WE ARE’ GOING TO SCAN ALL ELEMENTS AS IN SHCONN, 

GET THE STRESSES PERFORM ROTATION IF NECESSARY 

EVEU: NUM GF ROAS ELTSs LLIOSELT NUMBER 

NPREV:s NOOE NUM. (FOR EPSA) OF LAST POINT OF ROW JUST. 
BELOW THE ELEMENT ROA K; LLROW:s NUM. OF ELTS IN EACH ROW 


NPREV=0 
IA ELTS +1(FOR NOOES) +2 ARTIFICIAL NOOES FOR NPREV 
NSHAPE=4 
LLELZ=JNNI-JNOI 
00 200 K=1,LLEL 


LLROW=IA(JNGI+K-1) 
NPREV=NPREV+LLROW*+3 


00 100 J=t,LLROW 
LLIO=IACJNOBEG eK -1) +J-1 
OSCLIJSSTRECO*(LLID|1) +1) 
OS(2)=STRE(9*(LLIO“-1) +2) 
08(3)=0. 

5) WILL 8E THE STRESSES IN LOCAL COORD. 
0$(4)20S(1) 
0$(5)=0S(2) 


VON MISES AT TOP ANO 30TTOM OF SHELL: 


VON 


TOXY=STRE(9*(LLIO~-1)+3) 
OSX=STRE(9*e(LLIO|1) +4) 
OSY=STRE(9*(LLIO|1) +5) 
OSX=0SX*6./THICK 
OSY=OSY*6./THICK 
08xX1=0S(1)+0SxX 
OSY1=0S(2)+0SY 
0SxX2=0S(1)-0SX 
OSY2=0S(2) -0SY 


MISES CALCULATIONS MOOIFTIEO AFTERWARDS FOR MEMB.+BENOING STRESS 
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C NO SHEAR AT TOP OF SHELL,GET STRESSES AT TOP ANO BOTTOM ANO TAKE 
C MAX. VALUE 
GC 
c TOP: 
SIGM1=0SX1 
SIGM2=0SY1 
VONML=SORTC(SIGMI *SIGM1-SIGM1 *SIGM2+SIGM2*SIGM2) 
C SO0TTOM: 
SIGM1=0Sx2 
SIGM2=0SY2 


VONM2=SQRTCSIGM1 *SIGML<SIGMI*SIGM2+SIGM2sSIGM2) 


Cc 
VONM=AMAX1(VONM1, VONM?) 
Cc O.X. FOR VON “MISES (VONM) 
C IS THE SHELL CURVEO, NLELI'S ARE NOOES SURROUNOING ELT. (CEPSA) 
c 
TF CCACJVELO-1) NEO.) OR. (CAC JVELO=2) .NE.0.)) THEN 
NLELISJ+NPREV+1 
NLEL2=J +NPREV+2 
NLEL3SJ+1+NPREV+LLROWe 1 +2 
NLELUYSJ +NPREV+LLROWe 1 +2 
XSCACIXLOC+2*NLEL1-2) +A(CJXLOC+2eNLEL2<2))/u, 
XSX+(CACJXLOC +2eNLEL3—<-2) +AC JXLOC+2*eNLEL Une) )/4, 
YE(ACJXLOC+2eNLELlo1l)+AC JXLOC+2*eNLEL2-1))/4. 
YoY +(ACJXLOC +2 xNLEL3e1) +ACJXLOC+exNLELUnl))/4. 
CALL CHELEM(X,Y,0S) 
ENOTF 
Cc ° 
Cc REAOY TO OISPLAY RECORO 
c 
NRITE(LLS) LLIO,NSHAPE, (OISP(I),I=1,10), 
i CSTRACI), 1=1, 190), 0STREC1)-0S(5),0S8S(6), 
2 DSTRE(4),(OSCI),IL=1,3), 
3 (OSTREC(I),1=8,10),VONM 
G WRITE(15,80) LLIO-NSHAPE, (OS(I),-I21,5),VONM 
100 CONTINUE 


200 CONTINUE 
CALL CLOSEC(LLS) 
80 FORMAT(214,6E15.8) 


RETURN 

ENO 

SUBROUTINE CHELEM(X,Y,05) REPORT 

OIMENSION A(1) BLANC 

COMMON [AC1) BLANC 

EQUIVALENCE (CTAC1), ACL) ) BLANC 

COMMON /CPARA/ NSTEPS , NOEG , NENO , NSHEET , N2BO , N3B01 , CPARA 

1 N3802, INTRVL, OELT, NHTOT, NJOIN, NRELAX, ALPHA CPARA 

$ LEN S&8x MSPARA 

COMMON / SSIZE / I8GC1),-NQJ,NELTOT,NI180,NLOAO,NBRECT, SSIZE 

i NBQUAO, ISHEET, NPRPTS, NSPTS, NSTRPTS, NVPTS, NHPTS, SoG 

2 NSSTYP, NNJ, NNTOT, LGOSP, LIQUO, LBCALC SS LZE 
C REPORT 

COMMON /STAB / I8T(1)-JSSIZE,JSPAR,JVELO, JSTRE,JSJXMAS,JIELM,JOMAT, STAB 

1 JL180,JLOOP,JPRET,JLHIS,JSTRN,-JFORC, JXLOC, JNGI,JNNI,JNGSEG, STAB 

$ JLSTOE,JIELMCL, JSTIF,JOEFL,JFORLG, STAB 

2 JIFPAR,JFLPAR,JXCOORO,JYCOORO, JNELTAX, JOELTAY,JVMA,JSEFX, STAB 

3 JFLUFR, JPRINC, JVELRAO, JGENFR, JPRES, JCSEP STAB 


OIMENSION XMAT(3,3) 
DIMENSION OS(3) 
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CU EWNODEWNMWW WY EW 





00 tO [=1,3 

00 10 JjJ=1,3 

XMATCI,J)=0. 
10 CONTINUE 


inti se ROUTINE SE HANGES THE STRESSE IN AN ELEMENT INTO A 
IN A GLOBAL SECTANGULAR SYSTEM ,GIVEN LOCAL COORO. OF 
CENTROIO X,Y 


le CURVATURE IN Y OIR. IS NOMSZERO ,CALCULATE 
THE ROTATION MATRIX AT EACH POINT 
fe > CA VECO-1)SNE.0.) THEN 
RCOURY=1 ./4(JVELO=1)} 
THETAL=Y/RCOURY 
XMAT(1I,1)=1. 
XMAT(2,2)=COS(THETAL) 
XMAT(2,3)=SINC(THETA1L) 
XMAT( 3,3) =SXMAT (2-2) 
XMAT( 3,2) ="1.*XMAT (2,35) 
CALL PROOC(XMAT,OS(1)-05(2),0S(3)) 
ENOIF , 


AANAAARAANNA 


c 
C IF CURVATURE IN X OIRECTION IS NON ZERO: 
IF CACIVELO=2)-NE.0.) THEN 
OO" 20 [=1,3 
00 20 J=1,3 
XMAT(I..J)=0. 
20 CONTINUE 
RCOURX=1./A(JVELO=2) 
THETA2=xX/RCOURX 
XMAT(2,2)21. 
XMAT(1,1)=2COS(THETA2) 
XMAT(3,3)2XMAT(1,1) 
XMAT( 3,1) Se1.*2SIN(THETA2]) 
XMAT(L-3)2SINCTHETA2?) 
CALL PROOD(XMAT,0S(1),05(2),08S(3)) 
ENOIF 
RETURN 
ENO 
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